source: stable/1.17/Clp/src/unitTest.cpp

Last change on this file was 2531, checked in by stefan, 2 months ago

merge r2503,r2504,r2506,r2509..r2512,r2527 from trunk

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 93.2 KB
Line 
1/* $Id: unitTest.cpp 2531 2019-10-03 13:48:16Z 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#ifdef NDEBUG
7#undef NDEBUG
8#endif
9
10#include "ClpConfig.h"
11#include "CoinPragma.hpp"
12#include <cassert>
13#include <cstdio>
14#include <cstdlib>
15#include <cmath>
16#include <cfloat>
17#include <string>
18#include <iostream>
19
20#include "CoinMpsIO.hpp"
21#include "CoinPackedMatrix.hpp"
22#include "CoinPackedVector.hpp"
23#include "CoinStructuredModel.hpp"
24#include "CoinHelperFunctions.hpp"
25#include "CoinTime.hpp"
26#include "CoinFloatEqual.hpp"
27
28#if CLP_HAS_ABC
29#include "CoinAbcCommon.hpp"
30#endif
31#ifdef ABC_INHERIT
32#include "CoinAbcFactorization.hpp"
33#endif
34#include "ClpFactorization.hpp"
35#include "ClpSimplex.hpp"
36#include "ClpSimplexOther.hpp"
37#include "ClpSimplexNonlinear.hpp"
38#include "ClpInterior.hpp"
39#include "ClpLinearObjective.hpp"
40#include "ClpDualRowSteepest.hpp"
41#include "ClpDualRowDantzig.hpp"
42#include "ClpPrimalColumnSteepest.hpp"
43#include "ClpPrimalColumnDantzig.hpp"
44#include "ClpParameters.hpp"
45#include "ClpNetworkMatrix.hpp"
46#include "ClpPlusMinusOneMatrix.hpp"
47#include "MyMessageHandler.hpp"
48#include "MyEventHandler.hpp"
49
50#include "ClpPresolve.hpp"
51#include "Idiot.hpp"
52#if FACTORIZATION_STATISTICS
53extern double ftranTwiddleFactor1X;
54extern double ftranTwiddleFactor2X;
55extern double ftranFTTwiddleFactor1X;
56extern double ftranFTTwiddleFactor2X;
57extern double btranTwiddleFactor1X;
58extern double btranTwiddleFactor2X;
59extern double ftranFullTwiddleFactor1X;
60extern double ftranFullTwiddleFactor2X;
61extern double btranFullTwiddleFactor1X;
62extern double btranFullTwiddleFactor2X;
63extern double denseThresholdX;
64extern double twoThresholdX;
65extern double minRowsSparse;
66extern double largeRowsSparse;
67extern double mediumRowsDivider;
68extern double mediumRowsMinCount;
69extern double largeRowsCount;
70#endif
71
72//#############################################################################
73
74// Function Prototypes. Function definitions is in this file.
75static void testingMessage(const char *const msg);
76#if defined(COIN_HAS_AMD) || defined(COIN_HAS_CHOLMOD) || defined(COIN_HAS_GLPK)
77static int barrierAvailable = 1;
78static std::string nameBarrier = "barrier-UFL";
79#elif COIN_HAS_WSMP
80static int barrierAvailable = 2;
81static std::string nameBarrier = "barrier-WSSMP";
82#elif defined(COIN_HAS_MUMPS)
83static int barrierAvailable = 3;
84static std::string nameBarrier = "barrier-MUMPS";
85#else
86static int barrierAvailable = 0;
87static std::string nameBarrier = "barrier-slow";
88#endif
89#define NUMBER_ALGORITHMS 12
90// If you just want a subset then set some to 1
91static int switchOff[NUMBER_ALGORITHMS] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
92// shortName - 0 no , 1 yes
93ClpSolve setupForSolve(int algorithm, std::string &nameAlgorithm,
94  int shortName)
95{
96  ClpSolve solveOptions;
97  /* algorithms are
98        0 barrier
99        1 dual with volumne crash
100        2,3 dual with and without crash
101        4,5 primal with and without
102        6,7 automatic with and without
103        8,9 primal with idiot 1 and 5
104        10,11 primal with 70, dual with volume
105     */
106  switch (algorithm) {
107  case 0:
108    if (shortName)
109      nameAlgorithm = "ba";
110    else
111      nameAlgorithm = "nameBarrier";
112    solveOptions.setSolveType(ClpSolve::useBarrier);
113    if (barrierAvailable == 1)
114      solveOptions.setSpecialOption(4, 4);
115    else if (barrierAvailable == 2)
116      solveOptions.setSpecialOption(4, 2);
117    break;
118  case 1:
119#ifdef COIN_HAS_VOL
120    if (shortName)
121      nameAlgorithm = "du-vol-50";
122    else
123      nameAlgorithm = "dual-volume-50";
124    solveOptions.setSolveType(ClpSolve::useDual);
125    solveOptions.setSpecialOption(0, 2, 50); // volume
126#else
127    solveOptions.setSolveType(ClpSolve::notImplemented);
128#endif
129    break;
130  case 2:
131    if (shortName)
132      nameAlgorithm = "du-cr";
133    else
134      nameAlgorithm = "dual-crash";
135    solveOptions.setSolveType(ClpSolve::useDual);
136    solveOptions.setSpecialOption(0, 1);
137    break;
138  case 3:
139    if (shortName)
140      nameAlgorithm = "du";
141    else
142      nameAlgorithm = "dual";
143    solveOptions.setSolveType(ClpSolve::useDual);
144    break;
145  case 4:
146    if (shortName)
147      nameAlgorithm = "pr-cr";
148    else
149      nameAlgorithm = "primal-crash";
150    solveOptions.setSolveType(ClpSolve::usePrimal);
151    solveOptions.setSpecialOption(1, 1);
152    break;
153  case 5:
154    if (shortName)
155      nameAlgorithm = "pr";
156    else
157      nameAlgorithm = "primal";
158    solveOptions.setSolveType(ClpSolve::usePrimal);
159    break;
160  case 6:
161    if (shortName)
162      nameAlgorithm = "au-cr";
163    else
164      nameAlgorithm = "either-crash";
165    solveOptions.setSolveType(ClpSolve::automatic);
166    solveOptions.setSpecialOption(1, 1);
167    break;
168  case 7:
169    if (shortName)
170      nameAlgorithm = "au";
171    else
172      nameAlgorithm = "either";
173    solveOptions.setSolveType(ClpSolve::automatic);
174    break;
175  case 8:
176    if (shortName)
177      nameAlgorithm = "pr-id-1";
178    else
179      nameAlgorithm = "primal-idiot-1";
180    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
181    solveOptions.setSpecialOption(1, 2, 1); // idiot
182    break;
183  case 9:
184    if (shortName)
185      nameAlgorithm = "pr-id-5";
186    else
187      nameAlgorithm = "primal-idiot-5";
188    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
189    solveOptions.setSpecialOption(1, 2, 5); // idiot
190    break;
191  case 10:
192    if (shortName)
193      nameAlgorithm = "pr-id-70";
194    else
195      nameAlgorithm = "primal-idiot-70";
196    solveOptions.setSolveType(ClpSolve::usePrimalorSprint);
197    solveOptions.setSpecialOption(1, 2, 70); // idiot
198    break;
199  case 11:
200#ifdef COIN_HAS_VOL
201    if (shortName)
202      nameAlgorithm = "du-vol";
203    else
204      nameAlgorithm = "dual-volume";
205    solveOptions.setSolveType(ClpSolve::useDual);
206    solveOptions.setSpecialOption(0, 2, 3000); // volume
207#else
208    solveOptions.setSolveType(ClpSolve::notImplemented);
209#endif
210    break;
211  default:
212    abort();
213  }
214  if (shortName) {
215    // can switch off
216    if (switchOff[algorithm])
217      solveOptions.setSolveType(ClpSolve::notImplemented);
218  }
219  return solveOptions;
220}
221static void printSol(ClpSimplex &model)
222{
223  int numberRows = model.numberRows();
224  int numberColumns = model.numberColumns();
225
226  double *rowPrimal = model.primalRowSolution();
227  double *rowDual = model.dualRowSolution();
228  double *rowLower = model.rowLower();
229  double *rowUpper = model.rowUpper();
230
231  int iRow;
232  double objValue = model.getObjValue();
233  printf("Objvalue %g Rows (%d)\n", objValue, numberRows);
234  for (iRow = 0; iRow < numberRows; iRow++) {
235    printf("%d primal %g dual %g low %g up %g\n",
236      iRow, rowPrimal[iRow], rowDual[iRow],
237      rowLower[iRow], rowUpper[iRow]);
238  }
239  double *columnPrimal = model.primalColumnSolution();
240  double *columnDual = model.dualColumnSolution();
241  double *columnLower = model.columnLower();
242  double *columnUpper = model.columnUpper();
243  double offset;
244  //const double * gradient = model.objectiveAsObject()->gradient(&model,
245  //                                                       columnPrimal,offset,true,1);
246  const double *gradient = model.objective(columnPrimal, offset);
247  int iColumn;
248  objValue = -offset - model.objectiveOffset();
249  printf("offset %g (%g)\n", offset, model.objectiveOffset());
250  printf("Columns (%d)\n", numberColumns);
251  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
252    printf("%d primal %g dual %g low %g up %g\n",
253      iColumn, columnPrimal[iColumn], columnDual[iColumn],
254      columnLower[iColumn], columnUpper[iColumn]);
255    objValue += columnPrimal[iColumn] * gradient[iColumn];
256    if (fabs(columnPrimal[iColumn] * gradient[iColumn]) > 1.0e-8)
257      printf("obj -> %g gradient %g\n", objValue, gradient[iColumn]);
258  }
259  printf("Computed objective %g\n", objValue);
260}
261
262void usage(const std::string &key)
263{
264  std::cerr
265    << "Undefined parameter \"" << key << "\".\n"
266    << "Correct usage: \n"
267    << "  clp [-dirSample=V1] [-dirNetlib=V2] [-netlib]\n"
268    << "  where:\n"
269    << "    -dirSample: directory containing mps test files\n"
270    << "        Default value V1=\"../../Data/Sample\"\n"
271    << "    -dirNetlib: directory containing netlib files\"\n"
272    << "        Default value V2=\"../../Data/Netlib\"\n"
273    << "    -netlib\n"
274    << "        If specified, then netlib testset run as well as the nitTest.\n";
275}
276#if FACTORIZATION_STATISTICS
277int loSizeX = -1;
278int hiSizeX = 1000000;
279#endif
280//----------------------------------------------------------------
281#ifndef ABC_INHERIT
282#define AnySimplex ClpSimplex
283#else
284#include "AbcSimplex.hpp"
285#define AnySimplex AbcSimplex
286#endif
287int mainTest(int argc, const char *argv[], int algorithm,
288  AnySimplex empty, ClpSolve solveOptionsIn,
289  int switchOffValue, bool doVector)
290{
291  int i;
292
293  if (switchOffValue > 0) {
294    // switch off some
295    int iTest;
296    for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
297#ifndef NDEBUG
298      int bottom = switchOffValue % 10;
299      assert(bottom == 0 || bottom == 1);
300#endif
301      switchOffValue /= 10;
302      switchOff[iTest] = 0;
303    }
304  }
305  int numberFailures = 0;
306  // define valid parameter keywords
307  std::set< std::string > definedKeyWords;
308  definedKeyWords.insert("-dirSample");
309  definedKeyWords.insert("-dirNetlib");
310  definedKeyWords.insert("-netlib");
311
312  // Create a map of parameter keys and associated data
313  std::map< std::string, std::string > parms;
314  for (i = 1; i < argc; i++) {
315    std::string parm(argv[i]);
316    std::string key, value;
317    std::string::size_type eqPos = parm.find('=');
318
319    // Does parm contain and '='
320    if (eqPos == std::string::npos) {
321      //Parm does not contain '='
322      key = parm;
323    } else {
324      key = parm.substr(0, eqPos);
325      value = parm.substr(eqPos + 1);
326    }
327
328    // Is specifed key valid?
329    if (definedKeyWords.find(key) == definedKeyWords.end()) {
330      // invalid key word.
331      // Write help text
332      usage(key);
333      return 1;
334    }
335    parms[key] = value;
336  }
337
338  const char dirsep = CoinFindDirSeparator();
339  // Set directory containing mps data files.
340  std::string dirSample;
341  if (parms.find("-dirSample") != parms.end())
342    dirSample = parms["-dirSample"];
343  else
344    dirSample = dirsep == '/' ? "../../Data/Sample/" : "..\\..\\Data\\Sample\\";
345
346  // Set directory containing netlib data files.
347  std::string dirNetlib;
348  if (parms.find("-dirNetlib") != parms.end())
349    dirNetlib = parms["-dirNetlib"];
350  else
351    dirNetlib = dirsep == '/' ? "../../Data/Netlib/" : "..\\..\\Data\\Netlib\\";
352#if FACTORIZATION_STATISTICS == 0
353  if (!empty.numberRows()) {
354    testingMessage("Testing ClpSimplex\n");
355    ClpSimplexUnitTest(dirSample);
356  }
357#endif
358  if (parms.find("-netlib") != parms.end() || empty.numberRows()) {
359    unsigned int m;
360    std::string sizeLoHi;
361#if FACTORIZATION_STATISTICS
362    double ftranTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
363    double ftranTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
364    double ftranFTTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
365    double ftranFTTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
366    double btranTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
367    double btranTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
368    double ftranFullTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
369    double ftranFullTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
370    double btranFullTwiddleFactor1XChoice[3] = { 0.9, 1.0, 1.1 };
371    double btranFullTwiddleFactor2XChoice[3] = { 0.9, 1.0, 1.1 };
372    double denseThresholdXChoice[3] = { 2, 20, 40 };
373    double twoThresholdXChoice[3] = { 800, 1000, 1200 };
374    double minRowsSparseChoice[3] = { 250, 300, 350 };
375    double largeRowsSparseChoice[3] = { 8000, 10000, 12000 };
376    double mediumRowsDividerChoice[3] = { 5, 6, 7 };
377    double mediumRowsMinCountChoice[3] = { 300, 500, 600 };
378    double largeRowsCountChoice[3] = { 700, 1000, 1300 };
379    double times[3 * 3 * 3 * 3 * 3];
380    memset(times, 0, sizeof(times));
381#define whichParam(za, zname)            \
382  const char *choice##za = #zname;       \
383  const double *use##za = zname##Choice; \
384  double *external##za = &zname
385    whichParam(A, denseThresholdX);
386    whichParam(B, twoThresholdX);
387    whichParam(C, minRowsSparse);
388    whichParam(D, mediumRowsDivider);
389    whichParam(E, mediumRowsMinCount);
390#endif
391    // Define test problems:
392    //   mps names,
393    //   maximization or minimization,
394    //   Number of rows and columns in problem, and
395    //   objective function value
396    std::vector< std::string > mpsName;
397    std::vector< bool > min;
398    std::vector< int > nRows;
399    std::vector< int > nCols;
400    std::vector< double > objValue;
401    std::vector< double > objValueTol;
402    // 100 added means no presolve
403    std::vector< int > bestStrategy;
404    if (empty.numberRows()) {
405      std::string alg;
406      for (int iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
407        ClpSolve solveOptions = setupForSolve(iTest, alg, 0);
408        printf("%d %s ", iTest, alg.c_str());
409        if (switchOff[iTest])
410          printf("skipped by user\n");
411        else if (solveOptions.getSolveType() == ClpSolve::notImplemented)
412          printf("skipped as not available\n");
413        else
414          printf("will be tested\n");
415      }
416    }
417    if (!empty.numberRows()) {
418      mpsName.push_back("25fv47");
419      min.push_back(true);
420      nRows.push_back(822);
421      nCols.push_back(1571);
422      objValueTol.push_back(1.e-8);
423      objValue.push_back(5.5018458883E+03);
424      bestStrategy.push_back(0);
425      mpsName.push_back("80bau3b");
426      min.push_back(true);
427      nRows.push_back(2263);
428      nCols.push_back(9799);
429      objValueTol.push_back(1.e-8);
430      objValue.push_back(9.8722419241E+05);
431      bestStrategy.push_back(3);
432      mpsName.push_back("blend");
433      min.push_back(true);
434      nRows.push_back(75);
435      nCols.push_back(83);
436      objValueTol.push_back(1.e-8);
437      objValue.push_back(-3.0812149846e+01);
438      bestStrategy.push_back(3);
439      mpsName.push_back("pilotnov");
440      min.push_back(true);
441      nRows.push_back(976);
442      nCols.push_back(2172);
443      objValueTol.push_back(1.e-8);
444      objValue.push_back(-4.4972761882e+03);
445      bestStrategy.push_back(3);
446      mpsName.push_back("maros-r7");
447      min.push_back(true);
448      nRows.push_back(3137);
449      nCols.push_back(9408);
450      objValueTol.push_back(1.e-8);
451      objValue.push_back(1.4971851665e+06);
452      bestStrategy.push_back(2);
453      mpsName.push_back("pilot");
454      min.push_back(true);
455      nRows.push_back(1442);
456      nCols.push_back(3652);
457      objValueTol.push_back(1.e-5);
458      objValue.push_back(/*-5.5740430007e+02*/ -557.48972927292);
459      bestStrategy.push_back(3);
460      mpsName.push_back("pilot4");
461      min.push_back(true);
462      nRows.push_back(411);
463      nCols.push_back(1000);
464      objValueTol.push_back(5.e-5);
465      objValue.push_back(-2.5811392641e+03);
466      bestStrategy.push_back(3);
467      mpsName.push_back("pilot87");
468      min.push_back(true);
469      nRows.push_back(2031);
470      nCols.push_back(4883);
471      objValueTol.push_back(1.e-4);
472      objValue.push_back(3.0171072827e+02);
473      bestStrategy.push_back(0);
474      mpsName.push_back("adlittle");
475      min.push_back(true);
476      nRows.push_back(57);
477      nCols.push_back(97);
478      objValueTol.push_back(1.e-8);
479      objValue.push_back(2.2549496316e+05);
480      bestStrategy.push_back(3);
481      mpsName.push_back("afiro");
482      min.push_back(true);
483      nRows.push_back(28);
484      nCols.push_back(32);
485      objValueTol.push_back(1.e-8);
486      objValue.push_back(-4.6475314286e+02);
487      bestStrategy.push_back(3);
488      mpsName.push_back("agg");
489      min.push_back(true);
490      nRows.push_back(489);
491      nCols.push_back(163);
492      objValueTol.push_back(1.e-8);
493      objValue.push_back(-3.5991767287e+07);
494      bestStrategy.push_back(3);
495      mpsName.push_back("agg2");
496      min.push_back(true);
497      nRows.push_back(517);
498      nCols.push_back(302);
499      objValueTol.push_back(1.e-8);
500      objValue.push_back(-2.0239252356e+07);
501      bestStrategy.push_back(3);
502      mpsName.push_back("agg3");
503      min.push_back(true);
504      nRows.push_back(517);
505      nCols.push_back(302);
506      objValueTol.push_back(1.e-8);
507      objValue.push_back(1.0312115935e+07);
508      bestStrategy.push_back(4);
509      mpsName.push_back("bandm");
510      min.push_back(true);
511      nRows.push_back(306);
512      nCols.push_back(472);
513      objValueTol.push_back(1.e-8);
514      objValue.push_back(-1.5862801845e+02);
515      bestStrategy.push_back(2);
516      mpsName.push_back("beaconfd");
517      min.push_back(true);
518      nRows.push_back(174);
519      nCols.push_back(262);
520      objValueTol.push_back(1.e-8);
521      objValue.push_back(3.3592485807e+04);
522      bestStrategy.push_back(0);
523      mpsName.push_back("bnl1");
524      min.push_back(true);
525      nRows.push_back(644);
526      nCols.push_back(1175);
527      objValueTol.push_back(1.e-8);
528      objValue.push_back(1.9776295615E+03);
529      bestStrategy.push_back(3);
530      mpsName.push_back("bnl2");
531      min.push_back(true);
532      nRows.push_back(2325);
533      nCols.push_back(3489);
534      objValueTol.push_back(1.e-8);
535      objValue.push_back(1.8112365404e+03);
536      bestStrategy.push_back(3);
537      mpsName.push_back("boeing1");
538      min.push_back(true);
539      nRows.push_back(/*351*/ 352);
540      nCols.push_back(384);
541      objValueTol.push_back(1.e-8);
542      objValue.push_back(-3.3521356751e+02);
543      bestStrategy.push_back(3);
544      mpsName.push_back("boeing2");
545      min.push_back(true);
546      nRows.push_back(167);
547      nCols.push_back(143);
548      objValueTol.push_back(1.e-8);
549      objValue.push_back(-3.1501872802e+02);
550      bestStrategy.push_back(3);
551      mpsName.push_back("bore3d");
552      min.push_back(true);
553      nRows.push_back(234);
554      nCols.push_back(315);
555      objValueTol.push_back(1.e-8);
556      objValue.push_back(1.3730803942e+03);
557      bestStrategy.push_back(3);
558      mpsName.push_back("brandy");
559      min.push_back(true);
560      nRows.push_back(221);
561      nCols.push_back(249);
562      objValueTol.push_back(1.e-8);
563      objValue.push_back(1.5185098965e+03);
564      bestStrategy.push_back(3);
565      mpsName.push_back("capri");
566      min.push_back(true);
567      nRows.push_back(272);
568      nCols.push_back(353);
569      objValueTol.push_back(1.e-8);
570      objValue.push_back(2.6900129138e+03);
571      bestStrategy.push_back(3);
572      mpsName.push_back("cycle");
573      min.push_back(true);
574      nRows.push_back(1904);
575      nCols.push_back(2857);
576      objValueTol.push_back(1.e-9);
577      objValue.push_back(-5.2263930249e+00);
578      bestStrategy.push_back(3);
579      mpsName.push_back("czprob");
580      min.push_back(true);
581      nRows.push_back(930);
582      nCols.push_back(3523);
583      objValueTol.push_back(1.e-8);
584      objValue.push_back(2.1851966989e+06);
585      bestStrategy.push_back(3);
586      mpsName.push_back("d2q06c");
587      min.push_back(true);
588      nRows.push_back(2172);
589      nCols.push_back(5167);
590      objValueTol.push_back(1.e-7);
591      objValue.push_back(122784.21557456);
592      bestStrategy.push_back(0);
593      mpsName.push_back("d6cube");
594      min.push_back(true);
595      nRows.push_back(416);
596      nCols.push_back(6184);
597      objValueTol.push_back(1.e-7);
598      objValue.push_back(3.1549166667e+02);
599      bestStrategy.push_back(3);
600      mpsName.push_back("degen2");
601      min.push_back(true);
602      nRows.push_back(445);
603      nCols.push_back(534);
604      objValueTol.push_back(1.e-8);
605      objValue.push_back(-1.4351780000e+03);
606      bestStrategy.push_back(3);
607      mpsName.push_back("degen3");
608      min.push_back(true);
609      nRows.push_back(1504);
610      nCols.push_back(1818);
611      objValueTol.push_back(1.e-8);
612      objValue.push_back(-9.8729400000e+02);
613      bestStrategy.push_back(2);
614      mpsName.push_back("dfl001");
615      min.push_back(true);
616      nRows.push_back(6072);
617      nCols.push_back(12230);
618      objValueTol.push_back(1.e-5);
619      objValue.push_back(1.1266396047E+07);
620      bestStrategy.push_back(5);
621      mpsName.push_back("e226");
622      min.push_back(true);
623      nRows.push_back(224);
624      nCols.push_back(282);
625      objValueTol.push_back(1.e-8);
626      objValue.push_back(-1.8751929066e+01 + 7.113);
627      bestStrategy.push_back(3); // The correct answer includes -7.113 term. This is a constant in the objective function. See line 1683 of the mps file.
628      mpsName.push_back("etamacro");
629      min.push_back(true);
630      nRows.push_back(401);
631      nCols.push_back(688);
632      objValueTol.push_back(1.e-6);
633      objValue.push_back(-7.5571521774e+02);
634      bestStrategy.push_back(3);
635      mpsName.push_back("fffff800");
636      min.push_back(true);
637      nRows.push_back(525);
638      nCols.push_back(854);
639      objValueTol.push_back(1.e-6);
640      objValue.push_back(5.5567961165e+05);
641      bestStrategy.push_back(3);
642      mpsName.push_back("finnis");
643      min.push_back(true);
644      nRows.push_back(498);
645      nCols.push_back(614);
646      objValueTol.push_back(1.e-6);
647      objValue.push_back(1.7279096547e+05);
648      bestStrategy.push_back(3);
649      mpsName.push_back("fit1d");
650      min.push_back(true);
651      nRows.push_back(25);
652      nCols.push_back(1026);
653      objValueTol.push_back(1.e-8);
654      objValue.push_back(-9.1463780924e+03);
655      bestStrategy.push_back(3 + 100);
656      mpsName.push_back("fit1p");
657      min.push_back(true);
658      nRows.push_back(628);
659      nCols.push_back(1677);
660      objValueTol.push_back(1.e-8);
661      objValue.push_back(9.1463780924e+03);
662      bestStrategy.push_back(5 + 100);
663      mpsName.push_back("fit2d");
664      min.push_back(true);
665      nRows.push_back(26);
666      nCols.push_back(10500);
667      objValueTol.push_back(1.e-8);
668      objValue.push_back(-6.8464293294e+04);
669      bestStrategy.push_back(3 + 100);
670      mpsName.push_back("fit2p");
671      min.push_back(true);
672      nRows.push_back(3001);
673      nCols.push_back(13525);
674      objValueTol.push_back(1.e-9);
675      objValue.push_back(6.8464293232e+04);
676      bestStrategy.push_back(5 + 100);
677      mpsName.push_back("forplan");
678      min.push_back(true);
679      nRows.push_back(162);
680      nCols.push_back(421);
681      objValueTol.push_back(1.e-6);
682      objValue.push_back(-6.6421873953e+02);
683      bestStrategy.push_back(3);
684      mpsName.push_back("ganges");
685      min.push_back(true);
686      nRows.push_back(1310);
687      nCols.push_back(1681);
688      objValueTol.push_back(1.e-5);
689      objValue.push_back(-1.0958636356e+05);
690      bestStrategy.push_back(3);
691      mpsName.push_back("gfrd-pnc");
692      min.push_back(true);
693      nRows.push_back(617);
694      nCols.push_back(1092);
695      objValueTol.push_back(1.e-8);
696      objValue.push_back(6.9022359995e+06);
697      bestStrategy.push_back(3);
698      mpsName.push_back("greenbea");
699      min.push_back(true);
700      nRows.push_back(2393);
701      nCols.push_back(5405);
702      objValueTol.push_back(1.e-8);
703      objValue.push_back(/*-7.2462405908e+07*/ -72555248.129846);
704      bestStrategy.push_back(3);
705      mpsName.push_back("greenbeb");
706      min.push_back(true);
707      nRows.push_back(2393);
708      nCols.push_back(5405);
709      objValueTol.push_back(1.e-8);
710      objValue.push_back(/*-4.3021476065e+06*/ -4302260.2612066);
711      bestStrategy.push_back(3);
712      mpsName.push_back("grow15");
713      min.push_back(true);
714      nRows.push_back(301);
715      nCols.push_back(645);
716      objValueTol.push_back(1.e-8);
717      objValue.push_back(-1.0687094129e+08);
718      bestStrategy.push_back(4 + 100);
719      mpsName.push_back("grow22");
720      min.push_back(true);
721      nRows.push_back(441);
722      nCols.push_back(946);
723      objValueTol.push_back(1.e-8);
724      objValue.push_back(-1.6083433648e+08);
725      bestStrategy.push_back(4 + 100);
726      mpsName.push_back("grow7");
727      min.push_back(true);
728      nRows.push_back(141);
729      nCols.push_back(301);
730      objValueTol.push_back(1.e-8);
731      objValue.push_back(-4.7787811815e+07);
732      bestStrategy.push_back(4 + 100);
733      mpsName.push_back("israel");
734      min.push_back(true);
735      nRows.push_back(175);
736      nCols.push_back(142);
737      objValueTol.push_back(1.e-8);
738      objValue.push_back(-8.9664482186e+05);
739      bestStrategy.push_back(2);
740      mpsName.push_back("kb2");
741      min.push_back(true);
742      nRows.push_back(44);
743      nCols.push_back(41);
744      objValueTol.push_back(1.e-8);
745      objValue.push_back(-1.7499001299e+03);
746      bestStrategy.push_back(3);
747      mpsName.push_back("lotfi");
748      min.push_back(true);
749      nRows.push_back(154);
750      nCols.push_back(308);
751      objValueTol.push_back(1.e-8);
752      objValue.push_back(-2.5264706062e+01);
753      bestStrategy.push_back(3);
754      mpsName.push_back("maros");
755      min.push_back(true);
756      nRows.push_back(847);
757      nCols.push_back(1443);
758      objValueTol.push_back(1.e-8);
759      objValue.push_back(-5.8063743701e+04);
760      bestStrategy.push_back(3);
761      mpsName.push_back("modszk1");
762      min.push_back(true);
763      nRows.push_back(688);
764      nCols.push_back(1620);
765      objValueTol.push_back(1.e-8);
766      objValue.push_back(3.2061972906e+02);
767      bestStrategy.push_back(3);
768      mpsName.push_back("nesm");
769      min.push_back(true);
770      nRows.push_back(663);
771      nCols.push_back(2923);
772      objValueTol.push_back(1.e-5);
773      objValue.push_back(1.4076073035e+07);
774      bestStrategy.push_back(2);
775      mpsName.push_back("perold");
776      min.push_back(true);
777      nRows.push_back(626);
778      nCols.push_back(1376);
779      objValueTol.push_back(1.e-6);
780      objValue.push_back(-9.3807580773e+03);
781      bestStrategy.push_back(3);
782      //mpsName.push_back("qap12");min.push_back(true);nRows.push_back(3193);nCols.push_back(8856);objValueTol.push_back(1.e-6);objValue.push_back(5.2289435056e+02);bestStrategy.push_back(3);
783      //mpsName.push_back("qap15");min.push_back(true);nRows.push_back(6331);nCols.push_back(22275);objValueTol.push_back(1.e-8);objValue.push_back(1.0409940410e+03);bestStrategy.push_back(3);
784      mpsName.push_back("recipe");
785      min.push_back(true);
786      nRows.push_back(92);
787      nCols.push_back(180);
788      objValueTol.push_back(1.e-8);
789      objValue.push_back(-2.6661600000e+02);
790      bestStrategy.push_back(3);
791      mpsName.push_back("sc105");
792      min.push_back(true);
793      nRows.push_back(106);
794      nCols.push_back(103);
795      objValueTol.push_back(1.e-8);
796      objValue.push_back(-5.2202061212e+01);
797      bestStrategy.push_back(3);
798      mpsName.push_back("sc205");
799      min.push_back(true);
800      nRows.push_back(206);
801      nCols.push_back(203);
802      objValueTol.push_back(1.e-8);
803      objValue.push_back(-5.2202061212e+01);
804      bestStrategy.push_back(3);
805      mpsName.push_back("sc50a");
806      min.push_back(true);
807      nRows.push_back(51);
808      nCols.push_back(48);
809      objValueTol.push_back(1.e-8);
810      objValue.push_back(-6.4575077059e+01);
811      bestStrategy.push_back(3);
812      mpsName.push_back("sc50b");
813      min.push_back(true);
814      nRows.push_back(51);
815      nCols.push_back(48);
816      objValueTol.push_back(1.e-8);
817      objValue.push_back(-7.0000000000e+01);
818      bestStrategy.push_back(3);
819      mpsName.push_back("scagr25");
820      min.push_back(true);
821      nRows.push_back(472);
822      nCols.push_back(500);
823      objValueTol.push_back(1.e-8);
824      objValue.push_back(-1.4753433061e+07);
825      bestStrategy.push_back(3);
826      mpsName.push_back("scagr7");
827      min.push_back(true);
828      nRows.push_back(130);
829      nCols.push_back(140);
830      objValueTol.push_back(1.e-6);
831      objValue.push_back(-2.3313892548e+06);
832      bestStrategy.push_back(3);
833      mpsName.push_back("scfxm1");
834      min.push_back(true);
835      nRows.push_back(331);
836      nCols.push_back(457);
837      objValueTol.push_back(1.e-8);
838      objValue.push_back(1.8416759028e+04);
839      bestStrategy.push_back(3);
840      mpsName.push_back("scfxm2");
841      min.push_back(true);
842      nRows.push_back(661);
843      nCols.push_back(914);
844      objValueTol.push_back(1.e-8);
845      objValue.push_back(3.6660261565e+04);
846      bestStrategy.push_back(3);
847      mpsName.push_back("scfxm3");
848      min.push_back(true);
849      nRows.push_back(991);
850      nCols.push_back(1371);
851      objValueTol.push_back(1.e-8);
852      objValue.push_back(5.4901254550e+04);
853      bestStrategy.push_back(3);
854      mpsName.push_back("scorpion");
855      min.push_back(true);
856      nRows.push_back(389);
857      nCols.push_back(358);
858      objValueTol.push_back(1.e-8);
859      objValue.push_back(1.8781248227e+03);
860      bestStrategy.push_back(3);
861      mpsName.push_back("scrs8");
862      min.push_back(true);
863      nRows.push_back(491);
864      nCols.push_back(1169);
865      objValueTol.push_back(1.e-5);
866      objValue.push_back(9.0429998619e+02);
867      bestStrategy.push_back(2);
868      mpsName.push_back("scsd1");
869      min.push_back(true);
870      nRows.push_back(78);
871      nCols.push_back(760);
872      objValueTol.push_back(1.e-8);
873      objValue.push_back(8.6666666743e+00);
874      bestStrategy.push_back(3 + 100);
875      mpsName.push_back("scsd6");
876      min.push_back(true);
877      nRows.push_back(148);
878      nCols.push_back(1350);
879      objValueTol.push_back(1.e-8);
880      objValue.push_back(5.0500000078e+01);
881      bestStrategy.push_back(3 + 100);
882      mpsName.push_back("scsd8");
883      min.push_back(true);
884      nRows.push_back(398);
885      nCols.push_back(2750);
886      objValueTol.push_back(1.e-7);
887      objValue.push_back(9.0499999993e+02);
888      bestStrategy.push_back(1 + 100);
889      mpsName.push_back("sctap1");
890      min.push_back(true);
891      nRows.push_back(301);
892      nCols.push_back(480);
893      objValueTol.push_back(1.e-8);
894      objValue.push_back(1.4122500000e+03);
895      bestStrategy.push_back(3);
896      mpsName.push_back("sctap2");
897      min.push_back(true);
898      nRows.push_back(1091);
899      nCols.push_back(1880);
900      objValueTol.push_back(1.e-8);
901      objValue.push_back(1.7248071429e+03);
902      bestStrategy.push_back(3);
903      mpsName.push_back("sctap3");
904      min.push_back(true);
905      nRows.push_back(1481);
906      nCols.push_back(2480);
907      objValueTol.push_back(1.e-8);
908      objValue.push_back(1.4240000000e+03);
909      bestStrategy.push_back(3);
910      mpsName.push_back("seba");
911      min.push_back(true);
912      nRows.push_back(516);
913      nCols.push_back(1028);
914      objValueTol.push_back(1.e-8);
915      objValue.push_back(1.5711600000e+04);
916      bestStrategy.push_back(3);
917      mpsName.push_back("share1b");
918      min.push_back(true);
919      nRows.push_back(118);
920      nCols.push_back(225);
921      objValueTol.push_back(1.e-8);
922      objValue.push_back(-7.6589318579e+04);
923      bestStrategy.push_back(3);
924      mpsName.push_back("share2b");
925      min.push_back(true);
926      nRows.push_back(97);
927      nCols.push_back(79);
928      objValueTol.push_back(1.e-8);
929      objValue.push_back(-4.1573224074e+02);
930      bestStrategy.push_back(3);
931      mpsName.push_back("shell");
932      min.push_back(true);
933      nRows.push_back(537);
934      nCols.push_back(1775);
935      objValueTol.push_back(1.e-8);
936      objValue.push_back(1.2088253460e+09);
937      bestStrategy.push_back(3);
938      mpsName.push_back("ship04l");
939      min.push_back(true);
940      nRows.push_back(403);
941      nCols.push_back(2118);
942      objValueTol.push_back(1.e-8);
943      objValue.push_back(1.7933245380e+06);
944      bestStrategy.push_back(3);
945      mpsName.push_back("ship04s");
946      min.push_back(true);
947      nRows.push_back(403);
948      nCols.push_back(1458);
949      objValueTol.push_back(1.e-8);
950      objValue.push_back(1.7987147004e+06);
951      bestStrategy.push_back(3);
952      mpsName.push_back("ship08l");
953      min.push_back(true);
954      nRows.push_back(779);
955      nCols.push_back(4283);
956      objValueTol.push_back(1.e-8);
957      objValue.push_back(1.9090552114e+06);
958      bestStrategy.push_back(3);
959      mpsName.push_back("ship08s");
960      min.push_back(true);
961      nRows.push_back(779);
962      nCols.push_back(2387);
963      objValueTol.push_back(1.e-8);
964      objValue.push_back(1.9200982105e+06);
965      bestStrategy.push_back(2);
966      mpsName.push_back("ship12l");
967      min.push_back(true);
968      nRows.push_back(1152);
969      nCols.push_back(5427);
970      objValueTol.push_back(1.e-8);
971      objValue.push_back(1.4701879193e+06);
972      bestStrategy.push_back(3);
973      mpsName.push_back("ship12s");
974      min.push_back(true);
975      nRows.push_back(1152);
976      nCols.push_back(2763);
977      objValueTol.push_back(1.e-8);
978      objValue.push_back(1.4892361344e+06);
979      bestStrategy.push_back(2);
980      mpsName.push_back("sierra");
981      min.push_back(true);
982      nRows.push_back(1228);
983      nCols.push_back(2036);
984      objValueTol.push_back(1.e-8);
985      objValue.push_back(1.5394362184e+07);
986      bestStrategy.push_back(3);
987      mpsName.push_back("stair");
988      min.push_back(true);
989      nRows.push_back(357);
990      nCols.push_back(467);
991      objValueTol.push_back(1.e-8);
992      objValue.push_back(-2.5126695119e+02);
993      bestStrategy.push_back(3);
994      mpsName.push_back("standata");
995      min.push_back(true);
996      nRows.push_back(360);
997      nCols.push_back(1075);
998      objValueTol.push_back(1.e-8);
999      objValue.push_back(1.2576995000e+03);
1000      bestStrategy.push_back(3);
1001      //mpsName.push_back("standgub");min.push_back(true);nRows.push_back(362);nCols.push_back(1184);objValueTol.push_back(1.e-8);objValue.push_back(1257.6995); bestStrategy.push_back(3);
1002      mpsName.push_back("standmps");
1003      min.push_back(true);
1004      nRows.push_back(468);
1005      nCols.push_back(1075);
1006      objValueTol.push_back(1.e-8);
1007      objValue.push_back(1.4060175000E+03);
1008      bestStrategy.push_back(3);
1009      mpsName.push_back("stocfor1");
1010      min.push_back(true);
1011      nRows.push_back(118);
1012      nCols.push_back(111);
1013      objValueTol.push_back(1.e-8);
1014      objValue.push_back(-4.1131976219E+04);
1015      bestStrategy.push_back(3);
1016      mpsName.push_back("stocfor2");
1017      min.push_back(true);
1018      nRows.push_back(2158);
1019      nCols.push_back(2031);
1020      objValueTol.push_back(1.e-8);
1021      objValue.push_back(-3.9024408538e+04);
1022      bestStrategy.push_back(3);
1023      //mpsName.push_back("stocfor3");min.push_back(true);nRows.push_back(16676);nCols.push_back(15695);objValueTol.push_back(1.e-8);objValue.push_back(-3.9976661576e+04);bestStrategy.push_back(3);
1024      //mpsName.push_back("truss");min.push_back(true);nRows.push_back(1001);nCols.push_back(8806);objValueTol.push_back(1.e-8);objValue.push_back(4.5881584719e+05);bestStrategy.push_back(3);
1025      mpsName.push_back("tuff");
1026      min.push_back(true);
1027      nRows.push_back(334);
1028      nCols.push_back(587);
1029      objValueTol.push_back(1.e-8);
1030      objValue.push_back(2.9214776509e-01);
1031      bestStrategy.push_back(3);
1032      mpsName.push_back("vtpbase");
1033      min.push_back(true);
1034      nRows.push_back(199);
1035      nCols.push_back(203);
1036      objValueTol.push_back(1.e-8);
1037      objValue.push_back(1.2983146246e+05);
1038      bestStrategy.push_back(3);
1039      mpsName.push_back("wood1p");
1040      min.push_back(true);
1041      nRows.push_back(245);
1042      nCols.push_back(2594);
1043      objValueTol.push_back(5.e-5);
1044      objValue.push_back(1.4429024116e+00);
1045      bestStrategy.push_back(3);
1046      mpsName.push_back("woodw");
1047      min.push_back(true);
1048      nRows.push_back(1099);
1049      nCols.push_back(8405);
1050      objValueTol.push_back(1.e-8);
1051      objValue.push_back(1.3044763331E+00);
1052      bestStrategy.push_back(3);
1053    } else {
1054      // Just testing one
1055      mpsName.push_back(empty.problemName());
1056      min.push_back(true);
1057      nRows.push_back(-1);
1058      nCols.push_back(-1);
1059      objValueTol.push_back(1.e-8);
1060      objValue.push_back(0.0);
1061      bestStrategy.push_back(0);
1062      int iTest;
1063      std::string alg;
1064      for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
1065        ClpSolve solveOptions = setupForSolve(iTest, alg, 0);
1066        printf("%d %s ", iTest, alg.c_str());
1067        if (switchOff[iTest])
1068          printf("skipped by user\n");
1069        else if (solveOptions.getSolveType() == ClpSolve::notImplemented)
1070          printf("skipped as not available\n");
1071        else
1072          printf("will be tested\n");
1073      }
1074    }
1075
1076    double timeTaken = 0.0;
1077    if (!barrierAvailable)
1078      switchOff[0] = 1;
1079    // Loop once for each Mps File
1080    for (m = 0; m < mpsName.size(); m++) {
1081#if FACTORIZATION_STATISTICS
1082      if (nRows[m] < loSizeX || nRows[m] >= hiSizeX) {
1083        std::cerr << "  skipping mps file: " << mpsName[m]
1084                  << " as " << nRows[m]
1085                  << " (" << m + 1 << " out of " << mpsName.size() << ")" << std::endl;
1086        continue;
1087      }
1088#endif
1089      std::cerr << "  processing mps file: " << mpsName[m]
1090                << " (" << m + 1 << " out of " << mpsName.size() << ")" << std::endl;
1091      AnySimplex solutionBase = empty;
1092      std::string fn = dirNetlib + mpsName[m];
1093      if (!empty.numberRows() || algorithm < 6) {
1094        // Read data mps file,
1095        CoinMpsIO mps;
1096        int nerrors = mps.readMps(fn.c_str(), "mps");
1097        if (nerrors) {
1098          std::cerr << "Error " << nerrors << " when reading model from "
1099                    << fn.c_str() << "! Aborting tests.\n";
1100          return 1;
1101        }
1102        solutionBase.loadProblem(*mps.getMatrixByCol(), mps.getColLower(),
1103          mps.getColUpper(),
1104          mps.getObjCoefficients(),
1105          mps.getRowLower(), mps.getRowUpper());
1106
1107        solutionBase.setDblParam(ClpObjOffset, mps.objectiveOffset());
1108      }
1109
1110      // Runs through strategies
1111      if (algorithm == 6 || algorithm == 7) {
1112        // algorithms tested are at top of file
1113        double testTime[NUMBER_ALGORITHMS];
1114        std::string alg[NUMBER_ALGORITHMS];
1115        int iTest;
1116        for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
1117          ClpSolve solveOptions = setupForSolve(iTest, alg[iTest], 1);
1118          if (solveOptions.getSolveType() != ClpSolve::notImplemented) {
1119            double time1 = CoinCpuTime();
1120            ClpSimplex solution = solutionBase;
1121            if (solution.maximumSeconds() < 0.0)
1122              solution.setMaximumSeconds(120.0);
1123            if (doVector) {
1124              ClpMatrixBase *matrix = solution.clpMatrix();
1125              if (dynamic_cast< ClpPackedMatrix * >(matrix)) {
1126                ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
1127                clpMatrix->makeSpecialColumnCopy();
1128              }
1129            }
1130            solution.initialSolve(solveOptions);
1131            double time2 = CoinCpuTime() - time1;
1132            testTime[iTest] = time2;
1133            printf("Finished %s Took %g seconds (%d iterations) - status %d\n",
1134              mpsName[m].c_str(), time2, solution.problemStatus(), solution.numberIterations());
1135            if (solution.problemStatus())
1136              testTime[iTest] = 1.0e20;
1137          } else {
1138            testTime[iTest] = 1.0e30;
1139          }
1140        }
1141        int iBest = -1;
1142        double dBest = 1.0e10;
1143        printf("%s", fn.c_str());
1144        for (iTest = 0; iTest < NUMBER_ALGORITHMS; iTest++) {
1145          if (testTime[iTest] < 1.0e30) {
1146            printf(" %s %g",
1147              alg[iTest].c_str(), testTime[iTest]);
1148            if (testTime[iTest] < dBest) {
1149              dBest = testTime[iTest];
1150              iBest = iTest;
1151            }
1152          }
1153        }
1154        printf("\n");
1155        if (iBest >= 0)
1156          printf("Best strategy for %s is %s (%d) which takes %g seconds\n",
1157            fn.c_str(), alg[iBest].c_str(), iBest, testTime[iBest]);
1158        else
1159          printf("No strategy finished in time limit\n");
1160        continue;
1161      }
1162      double time1 = CoinCpuTime();
1163      AnySimplex solution = solutionBase;
1164#if 0
1165               solution.setOptimizationDirection(-1);
1166               {
1167                    int j;
1168                    double * obj = solution.objective();
1169                    int n = solution.numberColumns();
1170                    for (j = 0; j < n; j++)
1171                         obj[j] *= -1.0;
1172               }
1173#endif
1174      ClpSolve::SolveType method;
1175      ClpSolve solveOptions = solveOptionsIn;
1176      std::string nameAlgorithm;
1177      if (algorithm != 5) {
1178        if (algorithm == 0) {
1179          method = ClpSolve::useDual;
1180          nameAlgorithm = "dual";
1181        } else if (algorithm == 1) {
1182          method = ClpSolve::usePrimalorSprint;
1183          nameAlgorithm = "primal";
1184        } else if (algorithm == 3) {
1185          method = ClpSolve::automatic;
1186          nameAlgorithm = "either";
1187        } else {
1188          nameAlgorithm = "barrier-slow";
1189#if defined(COIN_HAS_AMD) || defined(COIN_HAS_CHOLMOD) || defined(COIN_HAS_GLPK)
1190          solveOptions.setSpecialOption(4, 4);
1191          nameAlgorithm = "barrier-UFL";
1192#endif
1193#ifdef COIN_HAS_WSMP
1194          solveOptions.setSpecialOption(4, 2);
1195          nameAlgorithm = "barrier-WSSMP";
1196#endif
1197#ifdef COIN_HAS_MUMPS
1198          solveOptions.setSpecialOption(4, 6);
1199          nameAlgorithm = "barrier-MUMPS";
1200#endif
1201          method = ClpSolve::useBarrier;
1202        }
1203        solveOptions.setSolveType(method);
1204      } else {
1205        int iAlg = bestStrategy[m];
1206        int presolveOff = iAlg / 100;
1207        iAlg = iAlg % 100;
1208        if (!barrierAvailable && iAlg == 0) {
1209          if (nRows[m] != 2172)
1210            iAlg = 5; // try primal
1211          else
1212            iAlg = 3; // d2q06c
1213        }
1214        solveOptions = setupForSolve(iAlg, nameAlgorithm, 0);
1215        if (presolveOff)
1216          solveOptions.setPresolveType(ClpSolve::presolveOff);
1217      }
1218      if (doVector) {
1219        ClpMatrixBase *matrix = solution.clpMatrix();
1220        if (dynamic_cast< ClpPackedMatrix * >(matrix)) {
1221          ClpPackedMatrix *clpMatrix = dynamic_cast< ClpPackedMatrix * >(matrix);
1222          clpMatrix->makeSpecialColumnCopy();
1223        }
1224      }
1225#if FACTORIZATION_STATISTICS
1226      double timesOne[3 * 3 * 3 * 3 * 3];
1227      memset(timesOne, 0, sizeof(timesOne));
1228      int iterationsOne[3 * 3 * 3 * 3 * 3];
1229      memset(iterationsOne, 0, sizeof(iterationsOne));
1230      AnySimplex saveModel(solution);
1231      double time2;
1232#if 0
1233               solution.initialSolve(solveOptions);
1234               time2 = CoinCpuTime() - time1;
1235               timeTaken += time2;
1236               printf("%s took %g seconds using algorithm %s\n", fn.c_str(), time2, nameAlgorithm.c_str());
1237#endif
1238#define loA 1
1239#define loB 0
1240#define loC 0
1241#define loD 1
1242#define loE 1
1243#define hiA 2
1244#define hiB 1
1245#define hiC 1
1246#define hiD 2
1247#define hiE 2
1248      time2 = CoinCpuTime();
1249      for (int iA = loA; iA < hiA; iA++) {
1250        *externalA = useA[iA];
1251        for (int iB = loB; iB < hiB; iB++) {
1252          *externalB = useB[iB];
1253          for (int iC = loC; iC < hiC; iC++) {
1254            *externalC = useC[iC];
1255            for (int iD = loD; iD < hiD; iD++) {
1256              *externalD = useD[iD];
1257              for (int iE = loE; iE < hiE; iE++) {
1258                *externalE = useE[iE];
1259                solution = saveModel;
1260                solution.initialSolve(solveOptions);
1261                double time3 = CoinCpuTime();
1262                printf("Finished %s Took %g seconds (%d iterations) - status %d\n",
1263                  mpsName[m].c_str(), time3 - time2, solution.numberIterations(), solution.problemStatus());
1264                iterationsOne[iA + 3 * iB + 9 * iC + 27 * iD + 81 * iE] = solution.numberIterations();
1265                timesOne[iA + 3 * iB + 9 * iC + 27 * iD + 81 * iE] = time3 - time2;
1266                times[iA + 3 * iB + 9 * iC + 27 * iD + 81 * iE] += time3 - time2;
1267                time2 = time3;
1268              }
1269            }
1270          }
1271        }
1272      }
1273      double bestTime = 1.0e100;
1274      int iBestTime = -1;
1275      double bestTimePer = 1.0e100;
1276      int iBestTimePer = -1;
1277      int bestIterations = 1000000;
1278      int iBestIterations = -1;
1279      printf("TTimes ");
1280      for (int i = 0; i < 3 * 3 * 3 * 3 * 3; i++) {
1281        // skip if not done
1282        if (!iterationsOne[i])
1283          continue;
1284        double average = timesOne[i] / static_cast< double >(iterationsOne[i]);
1285        if (timesOne[i] < bestTime) {
1286          bestTime = timesOne[i];
1287          iBestTime = i;
1288        }
1289        if (average < bestTimePer) {
1290          bestTimePer = average;
1291          iBestTimePer = i;
1292        }
1293        if (iterationsOne[i] < bestIterations) {
1294          bestIterations = iterationsOne[i];
1295          iBestIterations = i;
1296        }
1297        printf("%.2f ", timesOne[i]);
1298      }
1299      printf("\n");
1300      int iA, iB, iC, iD, iE;
1301      iA = iBestTime;
1302      iE = iA / 81;
1303      iA -= 81 * iE;
1304      iD = iA / 27;
1305      iA -= 27 * iD;
1306      iC = iA / 9;
1307      iA -= 9 * iC;
1308      iB = iA / 3;
1309      iA -= 3 * iB;
1310      printf("Best time %.2f %s=%g %s=%g %s=%g %s=%g %s=%g\n", bestTime, choiceA, useA[iA],
1311        choiceB, useB[iB], choiceC, useC[iC], choiceD, useD[iD], choiceE, useE[iE]);
1312      iA = iBestTimePer;
1313      iE = iA / 81;
1314      iA -= 81 * iE;
1315      iD = iA / 27;
1316      iA -= 27 * iD;
1317      iC = iA / 9;
1318      iA -= 9 * iC;
1319      iB = iA / 3;
1320      iA -= 3 * iB;
1321      printf("Best time per iteration %g %s=%g %s=%g %s=%g %s=%g %s=%g\n", bestTimePer, choiceA, useA[iA],
1322        choiceB, useB[iB], choiceC, useC[iC], choiceD, useD[iD], choiceE, useE[iE]);
1323      iA = iBestIterations;
1324      iE = iA / 81;
1325      iA -= 81 * iE;
1326      iD = iA / 27;
1327      iA -= 27 * iD;
1328      iC = iA / 9;
1329      iA -= 9 * iC;
1330      iB = iA / 3;
1331      iA -= 3 * iB;
1332      printf("Best iterations %d %s=%g %s=%g %s=%g %s=%g %s=%g\n", bestIterations, choiceA, useA[iA],
1333        choiceB, useB[iB], choiceC, useC[iC], choiceD, useD[iD], choiceE, useE[iE]);
1334#else
1335      solution.initialSolve(solveOptions);
1336      double time2 = CoinCpuTime() - time1;
1337      timeTaken += time2;
1338      printf("%s took %g seconds using algorithm %s\n", fn.c_str(), time2, nameAlgorithm.c_str());
1339#endif
1340      // Test objective solution value
1341      {
1342        if (!solution.isProvenOptimal()) {
1343          std::cerr << "** NOT OPTIMAL ";
1344          numberFailures++;
1345        }
1346        double soln = solution.objectiveValue();
1347        CoinRelFltEq eq(objValueTol[m]);
1348        std::cerr << soln << ",  " << objValue[m] << " diff " << soln - objValue[m] << std::endl;
1349        if (!eq(soln, objValue[m])) {
1350          printf("** difference fails\n");
1351          numberFailures++;
1352        }
1353      }
1354    }
1355    printf("Total time %g seconds\n", timeTaken);
1356#if FACTORIZATION_STATISTICS
1357    double bestTime = 1.0e100;
1358    int iBestTime = -1;
1359    for (int i = 0; i < 3 * 3 * 3 * 3 * 3; i++) {
1360      if (times[i] && times[i] < bestTime) {
1361        bestTime = times[i];
1362        iBestTime = i;
1363      }
1364    }
1365    for (int iA = 0; iA < hiA; iA++) {
1366      for (int iB = 0; iB < hiB; iB++) {
1367        for (int iC = 0; iC < hiC; iC++) {
1368          for (int iD = loD; iD < hiD; iD++) {
1369            for (int iE = loE; iE < hiE; iE++) {
1370              int k = iA + 3 * iB + 9 * iC + 27 * iD + 81 * iE;
1371              double thisTime = times[k];
1372              if (thisTime) {
1373                if (k == iBestTime)
1374                  printf("** ");
1375                printf("%s=%g %s=%g %s=%g %s=%g %s=%g - time %.2f\n",
1376                  choiceA, useA[iA],
1377                  choiceB, useB[iB], choiceC, useC[iC],
1378                  choiceD, useD[iD], choiceE, useE[iE], thisTime);
1379              }
1380            }
1381          }
1382        }
1383      }
1384    }
1385    //exit(0);
1386#endif
1387  } else {
1388    testingMessage("***Skipped Testing on netlib    ***\n");
1389    testingMessage("***use -netlib to test class***\n");
1390  }
1391  if (!numberFailures)
1392    testingMessage("All tests completed successfully\n");
1393  else
1394    testingMessage("Some tests failed\n");
1395  return 0;
1396}
1397
1398// Display message on stdout and stderr
1399static void testingMessage(const char *const msg)
1400{
1401  std::cerr << msg;
1402  //cout <<endl <<"*****************************************"
1403  //     <<endl <<msg <<endl;
1404}
1405
1406//--------------------------------------------------------------------------
1407// test factorization methods and simplex method and simple barrier
1408void ClpSimplexUnitTest(const std::string &dirSample)
1409{
1410
1411  CoinRelFltEq eq(0.000001);
1412
1413  {
1414    ClpSimplex solution;
1415
1416    // matrix data
1417    //deliberate hiccup of 2 between 0 and 1
1418    CoinBigIndex start[5] = { 0, 4, 7, 8, 9 };
1419    int length[5] = { 2, 3, 1, 1, 1 };
1420    int rows[11] = { 0, 2, -1, -1, 0, 1, 2, 0, 1, 2 };
1421    double elements[11] = { 7.0, 2.0, 1.0e10, 1.0e10, -2.0, 1.0, -2.0, 1, 1, 1 };
1422    CoinPackedMatrix matrix(true, 3, 5, 8, elements, rows, start, length);
1423
1424    // rim data
1425    double objective[7] = { -4.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
1426    double rowLower[5] = { 14.0, 3.0, 3.0, 1.0e10, 1.0e10 };
1427    double rowUpper[5] = { 14.0, 3.0, 3.0, -1.0e10, -1.0e10 };
1428    double colLower[7] = { 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0 };
1429    double colUpper[7] = { 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, 100.0 };
1430
1431    // basis 1
1432    int rowBasis1[3] = { -1, -1, -1 };
1433    int colBasis1[5] = { 1, 1, -1, -1, 1 };
1434    solution.loadProblem(matrix, colLower, colUpper, objective,
1435      rowLower, rowUpper);
1436    int i;
1437    solution.createStatus();
1438    for (i = 0; i < 3; i++) {
1439      if (rowBasis1[i] < 0) {
1440        solution.setRowStatus(i, ClpSimplex::atLowerBound);
1441      } else {
1442        solution.setRowStatus(i, ClpSimplex::basic);
1443      }
1444    }
1445    for (i = 0; i < 5; i++) {
1446      if (colBasis1[i] < 0) {
1447        solution.setColumnStatus(i, ClpSimplex::atLowerBound);
1448      } else {
1449        solution.setColumnStatus(i, ClpSimplex::basic);
1450      }
1451    }
1452    solution.setLogLevel(3 + 4 + 8 + 16 + 32);
1453    solution.primal();
1454    for (i = 0; i < 3; i++) {
1455      if (rowBasis1[i] < 0) {
1456        solution.setRowStatus(i, ClpSimplex::atLowerBound);
1457      } else {
1458        solution.setRowStatus(i, ClpSimplex::basic);
1459      }
1460    }
1461    for (i = 0; i < 5; i++) {
1462      if (colBasis1[i] < 0) {
1463        solution.setColumnStatus(i, ClpSimplex::atLowerBound);
1464      } else {
1465        solution.setColumnStatus(i, ClpSimplex::basic);
1466      }
1467    }
1468    // intricate stuff does not work with scaling
1469    solution.scaling(0);
1470#ifndef NDEBUG
1471    int returnCode = solution.factorize();
1472    assert(!returnCode);
1473#else
1474    solution.factorize();
1475#endif
1476    const double *colsol = solution.primalColumnSolution();
1477    const double *rowsol = solution.primalRowSolution();
1478    solution.getSolution(rowsol, colsol);
1479#ifndef NDEBUG
1480    double colsol1[5] = { 20.0 / 7.0, 3.0, 0.0, 0.0, 23.0 / 7.0 };
1481    for (i = 0; i < 5; i++) {
1482      assert(eq(colsol[i], colsol1[i]));
1483    }
1484    // now feed in again without actually doing factorization
1485#endif
1486    ClpFactorization factorization2 = *solution.factorization();
1487    ClpSimplex solution2 = solution;
1488    solution2.setFactorization(factorization2);
1489    solution2.createStatus();
1490    for (i = 0; i < 3; i++) {
1491      if (rowBasis1[i] < 0) {
1492        solution2.setRowStatus(i, ClpSimplex::atLowerBound);
1493      } else {
1494        solution2.setRowStatus(i, ClpSimplex::basic);
1495      }
1496    }
1497    for (i = 0; i < 5; i++) {
1498      if (colBasis1[i] < 0) {
1499        solution2.setColumnStatus(i, ClpSimplex::atLowerBound);
1500      } else {
1501        solution2.setColumnStatus(i, ClpSimplex::basic);
1502      }
1503    }
1504    // intricate stuff does not work with scaling
1505    solution2.scaling(0);
1506    solution2.getSolution(rowsol, colsol);
1507    colsol = solution2.primalColumnSolution();
1508    rowsol = solution2.primalRowSolution();
1509    for (i = 0; i < 5; i++) {
1510      assert(eq(colsol[i], colsol1[i]));
1511    }
1512    solution2.setDualBound(0.1);
1513    solution2.dual();
1514    objective[2] = -1.0;
1515    objective[3] = -0.5;
1516    objective[4] = 10.0;
1517    solution.dual();
1518    for (i = 0; i < 3; i++) {
1519      rowLower[i] = -1.0e20;
1520      colUpper[i + 2] = 0.0;
1521    }
1522    solution.setLogLevel(3);
1523    solution.dual();
1524    double rowObjective[] = { 1.0, 0.5, -10.0 };
1525    solution.loadProblem(matrix, colLower, colUpper, objective,
1526      rowLower, rowUpper, rowObjective);
1527    solution.dual();
1528    solution.loadProblem(matrix, colLower, colUpper, objective,
1529      rowLower, rowUpper, rowObjective);
1530    solution.primal();
1531  }
1532#ifndef COIN_NO_CLP_MESSAGE
1533  {
1534    CoinMpsIO m;
1535    std::string fn = dirSample + "exmip1";
1536    if (m.readMps(fn.c_str(), "mps") == 0) {
1537      ClpSimplex solution;
1538      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
1539        m.getObjCoefficients(),
1540        m.getRowLower(), m.getRowUpper());
1541      solution.dual();
1542      // Test event handling
1543      MyEventHandler handler;
1544      solution.passInEventHandler(&handler);
1545      int numberRows = solution.numberRows();
1546      // make sure values pass has something to do
1547      for (int i = 0; i < numberRows; i++)
1548        solution.setRowStatus(i, ClpSimplex::basic);
1549      solution.objective()[1] = -2.0;
1550      solution.primal(1);
1551      assert(solution.secondaryStatus() == 102); // Came out at end of pass
1552    } else {
1553      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
1554    }
1555  }
1556  // Test Message handler
1557  {
1558    CoinMpsIO m;
1559    std::string fn = dirSample + "exmip1";
1560    //fn = "Test/subGams4";
1561    if (m.readMps(fn.c_str(), "mps") == 0) {
1562      ClpSimplex model;
1563      model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
1564        m.getObjCoefficients(),
1565        m.getRowLower(), m.getRowUpper());
1566      // Message handler
1567      MyMessageHandler messageHandler(&model);
1568      std::cout << "Testing derived message handler" << std::endl;
1569      model.passInMessageHandler(&messageHandler);
1570      model.messagesPointer()->setDetailMessage(1, 102);
1571      model.setFactorizationFrequency(10);
1572      model.primal();
1573      model.primal(0, 3);
1574      model.setObjCoeff(3, -2.9473684210526314);
1575      model.primal(0, 3);
1576      // Write saved solutions
1577      int nc = model.getNumCols();
1578      size_t s;
1579      std::deque< StdVectorDouble > fep = messageHandler.getFeasibleExtremePoints();
1580      size_t numSavedSolutions = fep.size();
1581      for (s = 0; s < numSavedSolutions; ++s) {
1582        const StdVectorDouble &solnVec = fep[s];
1583        for (int c = 0; c < nc; ++c) {
1584          if (fabs(solnVec[c]) > 1.0e-8)
1585            std::cout << "Saved Solution: " << s << " ColNum: " << c << " Value: " << solnVec[c] << std::endl;
1586        }
1587      }
1588      // Solve again without scaling
1589      // and maximize then minimize
1590      messageHandler.clearFeasibleExtremePoints();
1591      model.scaling(0);
1592      model.setOptimizationDirection(-1);
1593      model.primal();
1594      model.setOptimizationDirection(1);
1595      model.primal();
1596      fep = messageHandler.getFeasibleExtremePoints();
1597      numSavedSolutions = fep.size();
1598      for (s = 0; s < numSavedSolutions; ++s) {
1599        const StdVectorDouble &solnVec = fep[s];
1600        for (int c = 0; c < nc; ++c) {
1601          if (fabs(solnVec[c]) > 1.0e-8)
1602            std::cout << "Saved Solution: " << s << " ColNum: " << c << " Value: " << solnVec[c] << std::endl;
1603        }
1604      }
1605    } else {
1606      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
1607    }
1608  }
1609#endif
1610  // Test dual ranging
1611  {
1612    CoinMpsIO m;
1613    std::string fn = dirSample + "exmip1";
1614    if (m.readMps(fn.c_str(), "mps") == 0) {
1615      ClpSimplex model;
1616      model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
1617        m.getObjCoefficients(),
1618        m.getRowLower(), m.getRowUpper());
1619      model.primal();
1620      int which[13] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 };
1621      double costIncrease[13];
1622      int sequenceIncrease[13];
1623      double costDecrease[13];
1624      int sequenceDecrease[13];
1625      // ranging
1626      model.dualRanging(13, which, costIncrease, sequenceIncrease,
1627        costDecrease, sequenceDecrease);
1628      int i;
1629      for (i = 0; i < 13; i++)
1630        printf("%d increase %g %d, decrease %g %d\n",
1631          i, costIncrease[i], sequenceIncrease[i],
1632          costDecrease[i], sequenceDecrease[i]);
1633      assert(fabs(costDecrease[3]) < 1.0e-4);
1634      assert(fabs(costIncrease[7] - 1.0) < 1.0e-4);
1635      model.setOptimizationDirection(-1);
1636      {
1637        int j;
1638        double *obj = model.objective();
1639        int n = model.numberColumns();
1640        for (j = 0; j < n; j++)
1641          obj[j] *= -1.0;
1642      }
1643      double costIncrease2[13];
1644      int sequenceIncrease2[13];
1645      double costDecrease2[13];
1646      int sequenceDecrease2[13];
1647      // ranging
1648      model.dualRanging(13, which, costIncrease2, sequenceIncrease2,
1649        costDecrease2, sequenceDecrease2);
1650      for (i = 0; i < 13; i++) {
1651        assert(fabs(costIncrease[i] - costDecrease2[i]) < 1.0e-6);
1652        assert(fabs(costDecrease[i] - costIncrease2[i]) < 1.0e-6);
1653        assert(sequenceIncrease[i] == sequenceDecrease2[i]);
1654        assert(sequenceDecrease[i] == sequenceIncrease2[i]);
1655      }
1656      // Now delete all rows and see what happens
1657      model.deleteRows(model.numberRows(), which);
1658      model.primal();
1659      // ranging
1660      if (!model.dualRanging(8, which, costIncrease, sequenceIncrease,
1661            costDecrease, sequenceDecrease)) {
1662        for (i = 0; i < 8; i++) {
1663          printf("%d increase %g %d, decrease %g %d\n",
1664            i, costIncrease[i], sequenceIncrease[i],
1665            costDecrease[i], sequenceDecrease[i]);
1666        }
1667      }
1668    } else {
1669      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
1670    }
1671  }
1672  // Test primal ranging
1673  {
1674    CoinMpsIO m;
1675    std::string fn = dirSample + "exmip1";
1676    if (m.readMps(fn.c_str(), "mps") == 0) {
1677      ClpSimplex model;
1678      model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
1679        m.getObjCoefficients(),
1680        m.getRowLower(), m.getRowUpper());
1681      model.primal();
1682      int which[13] = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12 };
1683      double valueIncrease[13];
1684      int sequenceIncrease[13];
1685      double valueDecrease[13];
1686      int sequenceDecrease[13];
1687      // ranging
1688      model.primalRanging(13, which, valueIncrease, sequenceIncrease,
1689        valueDecrease, sequenceDecrease);
1690      int i;
1691      for (i = 0; i < 13; i++)
1692        printf("%d increase %g %d, decrease %g %d\n",
1693          i, valueIncrease[i], sequenceIncrease[i],
1694          valueDecrease[i], sequenceDecrease[i]);
1695      assert(fabs(valueIncrease[3] - 0.642857) < 1.0e-4);
1696      assert(fabs(valueIncrease[8] - 2.95113) < 1.0e-4);
1697    } else {
1698      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
1699    }
1700#if 0
1701          // out until I find optimization bug
1702          // Test parametrics
1703          ClpSimplexOther * model2 = (ClpSimplexOther *) (&model);
1704          double rhs[] = { 1.0, 2.0, 3.0, 4.0, 5.0};
1705          double endingTheta = 1.0;
1706          model2->scaling(0);
1707          model2->setLogLevel(63);
1708          model2->parametrics(0.0, endingTheta, 0.1,
1709                              NULL, NULL, rhs, rhs, NULL);
1710#endif
1711  }
1712  // Test binv etc
1713  {
1714    /*
1715             Wolsey : Page 130
1716             max 4x1 -  x2
1717             7x1 - 2x2    <= 14
1718             x2    <= 3
1719             2x1 - 2x2    <= 3
1720             x1 in Z+, x2 >= 0
1721
1722             note slacks are -1 in Clp so signs may be different
1723          */
1724
1725    int n_cols = 2;
1726    int n_rows = 3;
1727
1728    double obj[2] = { -4.0, 1.0 };
1729    double collb[2] = { 0.0, 0.0 };
1730    double colub[2] = { COIN_DBL_MAX, COIN_DBL_MAX };
1731    double rowlb[3] = { -COIN_DBL_MAX, -COIN_DBL_MAX, -COIN_DBL_MAX };
1732    double rowub[3] = { 14.0, 3.0, 3.0 };
1733
1734    int rowIndices[5] = { 0, 2, 0, 1, 2 };
1735    int colIndices[5] = { 0, 0, 1, 1, 1 };
1736    double elements[5] = { 7.0, 2.0, -2.0, 1.0, -2.0 };
1737    CoinPackedMatrix M(true, rowIndices, colIndices, elements, 5);
1738
1739    ClpSimplex model;
1740    model.loadProblem(M, collb, colub, obj, rowlb, rowub);
1741    model.dual(0, 1); // keep factorization
1742
1743    //check that the tableau matches wolsey (B-1 A)
1744    // slacks in second part of binvA
1745    double *binvA = reinterpret_cast< double * >(malloc((n_cols + n_rows) * sizeof(double)));
1746
1747    printf("B-1 A by row\n");
1748    int i;
1749    for (i = 0; i < n_rows; i++) {
1750      model.getBInvARow(i, binvA, binvA + n_cols);
1751      printf("row: %d -> ", i);
1752      for (int j = 0; j < n_cols + n_rows; j++) {
1753        printf("%g, ", binvA[j]);
1754      }
1755      printf("\n");
1756    }
1757    // See if can re-use factorization AND arrays
1758    model.primal(0, 3 + 4); // keep factorization
1759    // And do by column
1760    printf("B-1 A by column\n");
1761    for (i = 0; i < n_rows + n_cols; i++) {
1762      model.getBInvACol(i, binvA);
1763      printf("column: %d -> ", i);
1764      for (int j = 0; j < n_rows; j++) {
1765        printf("%g, ", binvA[j]);
1766      }
1767      printf("\n");
1768    }
1769    /* Do twice -
1770             without and with scaling
1771          */
1772    // set scaling off
1773    model.scaling(0);
1774    for (int iPass = 0; iPass < 2; iPass++) {
1775      model.primal(0, 3 + 4); // keep factorization
1776      const double *rowScale = model.rowScale();
1777      const double *columnScale = model.columnScale();
1778      if (!iPass)
1779        assert(!rowScale);
1780      else
1781        assert(rowScale); // only true for this example
1782      /* has to be exactly correct as in OsiClpsolverInterface.cpp
1783                  (also redo each pass as may change
1784               */
1785      printf("B-1 A");
1786      for (i = 0; i < n_rows; i++) {
1787        model.getBInvARow(i, binvA, binvA + n_cols);
1788        printf("\nrow: %d -> ", i);
1789        int j;
1790        // First columns
1791        for (j = 0; j < n_cols; j++) {
1792          if (binvA[j]) {
1793            printf("(%d %g), ", j, binvA[j]);
1794          }
1795        }
1796        // now rows
1797        for (j = 0; j < n_rows; j++) {
1798          if (binvA[j + n_cols]) {
1799            printf("(%d %g), ", j + n_cols, binvA[j + n_cols]);
1800          }
1801        }
1802      }
1803      printf("\n");
1804      printf("And by column (trickier)");
1805      const int *pivotVariable = model.pivotVariable();
1806      for (i = 0; i < n_cols + n_rows; i++) {
1807        model.getBInvACol(i, binvA);
1808        printf("\ncolumn: %d -> ", i);
1809        for (int j = 0; j < n_rows; j++) {
1810          if (binvA[j]) {
1811            // need to know pivot variable for +1/-1 (slack) and row/column scaling
1812            int pivot = pivotVariable[j];
1813            if (pivot < n_cols) {
1814              // scaled coding is in just in case
1815              if (!columnScale) {
1816                printf("(%d %g), ", j, binvA[j]);
1817              } else {
1818                printf("(%d %g), ", j, binvA[j] * columnScale[pivot]);
1819              }
1820            } else {
1821              if (!rowScale) {
1822                printf("(%d %g), ", j, binvA[j]);
1823              } else {
1824                printf("(%d %g), ", j, binvA[j] / rowScale[pivot - n_cols]);
1825              }
1826            }
1827          }
1828        }
1829      }
1830      printf("\n");
1831      printf("binvrow");
1832      for (i = 0; i < n_rows; i++) {
1833        model.getBInvRow(i, binvA);
1834        printf("\nrow: %d -> ", i);
1835        int j;
1836        for (j = 0; j < n_rows; j++) {
1837          if (binvA[j])
1838            printf("(%d %g), ", j, binvA[j]);
1839        }
1840      }
1841      printf("\n");
1842      printf("And by column ");
1843      for (i = 0; i < n_rows; i++) {
1844        model.getBInvCol(i, binvA);
1845        printf("\ncol: %d -> ", i);
1846        int j;
1847        for (j = 0; j < n_rows; j++) {
1848          if (binvA[j])
1849            printf("(%d %g), ", j, binvA[j]);
1850        }
1851      }
1852      printf("\n");
1853      // now deal with next pass
1854      if (!iPass) {
1855        // get scaling for testing
1856        model.scaling(1);
1857      }
1858    }
1859    free(binvA);
1860    model.setColUpper(1, 2.0);
1861    model.dual(0, 2 + 4); // use factorization and arrays
1862    model.dual(0, 2); // hopefully will not use factorization
1863    model.primal(0, 3 + 4); // keep factorization
1864    // but say basis has changed
1865    model.setWhatsChanged(model.whatsChanged() & (~512));
1866    model.dual(0, 2); // hopefully will not use factorization
1867  }
1868  // test steepest edge
1869  {
1870    CoinMpsIO m;
1871    std::string fn = dirSample + "finnis";
1872    int returnCode = m.readMps(fn.c_str(), "mps");
1873    if (returnCode) {
1874      // probable cause is that gz not there
1875      fprintf(stderr, "Unable to open finnis.mps in %s\n", dirSample.c_str());
1876      fprintf(stderr, "Most probable cause is that sample data is not available, or finnis.mps is gzipped i.e. finnis.mps.gz and libz has not been activated\n");
1877      fprintf(stderr, "Either gunzip files or edit Makefiles/Makefile.location to get libz\n");
1878    } else {
1879      ClpModel model;
1880      model.loadProblem(*m.getMatrixByCol(), m.getColLower(),
1881        m.getColUpper(),
1882        m.getObjCoefficients(),
1883        m.getRowLower(), m.getRowUpper());
1884      ClpSimplex solution(model);
1885
1886      solution.scaling(1);
1887      solution.setDualBound(1.0e8);
1888      //solution.factorization()->maximumPivots(1);
1889      //solution.setLogLevel(3);
1890      solution.setDualTolerance(1.0e-7);
1891      // set objective sense,
1892      ClpDualRowSteepest steep;
1893      solution.setDualRowPivotAlgorithm(steep);
1894      solution.setDblParam(ClpObjOffset, m.objectiveOffset());
1895      solution.dual();
1896    }
1897  }
1898  // test normal solution
1899  {
1900    CoinMpsIO m;
1901    std::string fn = dirSample + "afiro";
1902    if (m.readMps(fn.c_str(), "mps") == 0) {
1903      ClpSimplex solution;
1904      ClpModel model;
1905      // do twice - without and with scaling
1906      int iPass;
1907      for (iPass = 0; iPass < 2; iPass++) {
1908        // explicit row objective for testing
1909        int nr = m.getNumRows();
1910        double *rowObj = new double[nr];
1911        CoinFillN(rowObj, nr, 0.0);
1912        model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
1913          m.getObjCoefficients(),
1914          m.getRowLower(), m.getRowUpper(), rowObj);
1915        delete[] rowObj;
1916        solution = ClpSimplex(model);
1917        if (iPass) {
1918          solution.scaling();
1919        }
1920        solution.dual();
1921        solution.dual();
1922        // test optimal
1923        assert(solution.status() == 0);
1924        int numberColumns = solution.numberColumns();
1925        int numberRows = solution.numberRows();
1926        CoinPackedVector colsol(numberColumns, solution.primalColumnSolution());
1927        double *objective = solution.objective();
1928#ifndef NDEBUG
1929        double objValue = colsol.dotProduct(objective);
1930#endif
1931        CoinRelFltEq eq(1.0e-8);
1932        assert(eq(objValue, -4.6475314286e+02));
1933        solution.dual();
1934        assert(eq(solution.objectiveValue(), -4.6475314286e+02));
1935        double *lower = solution.columnLower();
1936        double *upper = solution.columnUpper();
1937        double *sol = solution.primalColumnSolution();
1938        double *result = new double[numberColumns];
1939        CoinFillN(result, numberColumns, 0.0);
1940        solution.matrix()->transposeTimes(solution.dualRowSolution(), result);
1941        int iRow, iColumn;
1942        // see if feasible and dual feasible
1943        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1944          double value = sol[iColumn];
1945          assert(value < upper[iColumn] + 1.0e-8);
1946          assert(value > lower[iColumn] - 1.0e-8);
1947          value = objective[iColumn] - result[iColumn];
1948          assert(value > -1.0e-5);
1949          if (sol[iColumn] > 1.0e-5)
1950            assert(value < 1.0e-5);
1951        }
1952        delete[] result;
1953        result = new double[numberRows];
1954        CoinFillN(result, numberRows, 0.0);
1955        solution.matrix()->times(colsol, result);
1956        lower = solution.rowLower();
1957        upper = solution.rowUpper();
1958        sol = solution.primalRowSolution();
1959#ifndef NDEBUG
1960        for (iRow = 0; iRow < numberRows; iRow++) {
1961          double value = result[iRow];
1962          assert(eq(value, sol[iRow]));
1963          assert(value < upper[iRow] + 1.0e-8);
1964          assert(value > lower[iRow] - 1.0e-8);
1965        }
1966#endif
1967        delete[] result;
1968        // test row objective
1969        double *rowObjective = solution.rowObjective();
1970        CoinDisjointCopyN(solution.dualRowSolution(), numberRows, rowObjective);
1971        CoinDisjointCopyN(solution.dualColumnSolution(), numberColumns, objective);
1972        // this sets up all slack basis
1973        solution.createStatus();
1974        solution.dual();
1975        CoinFillN(rowObjective, numberRows, 0.0);
1976        CoinDisjointCopyN(m.getObjCoefficients(), numberColumns, objective);
1977        solution.dual();
1978      }
1979    } else {
1980      std::cerr << "Error reading afiro from sample data. Skipping test." << std::endl;
1981    }
1982  }
1983  // test unbounded
1984  {
1985    CoinMpsIO m;
1986    std::string fn = dirSample + "brandy";
1987    if (m.readMps(fn.c_str(), "mps") == 0) {
1988      ClpSimplex solution;
1989      // do twice - without and with scaling
1990      int iPass;
1991      for (iPass = 0; iPass < 2; iPass++) {
1992        solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
1993          m.getObjCoefficients(),
1994          m.getRowLower(), m.getRowUpper());
1995        if (iPass)
1996          solution.scaling();
1997        solution.setOptimizationDirection(-1);
1998        // test unbounded and ray
1999#ifdef DUAL
2000        solution.setDualBound(100.0);
2001        solution.dual();
2002#else
2003        solution.primal();
2004#endif
2005        assert(solution.status() == 2);
2006        int numberColumns = solution.numberColumns();
2007        int numberRows = solution.numberRows();
2008        double *lower = solution.columnLower();
2009        double *upper = solution.columnUpper();
2010        double *sol = solution.primalColumnSolution();
2011        double *ray = solution.unboundedRay();
2012        double *objective = solution.objective();
2013        double objChange = 0.0;
2014        int iRow, iColumn;
2015        // make sure feasible and columns form ray
2016        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2017          double value = sol[iColumn];
2018          assert(value < upper[iColumn] + 1.0e-8);
2019          assert(value > lower[iColumn] - 1.0e-8);
2020          value = ray[iColumn];
2021          if (value > 0.0)
2022            assert(upper[iColumn] > 1.0e30);
2023          else if (value < 0.0)
2024            assert(lower[iColumn] < -1.0e30);
2025          objChange += value * objective[iColumn];
2026        }
2027        // make sure increasing objective
2028        assert(objChange > 0.0);
2029        double *result = new double[numberRows];
2030        CoinFillN(result, numberRows, 0.0);
2031        solution.matrix()->times(sol, result);
2032        lower = solution.rowLower();
2033        upper = solution.rowUpper();
2034        sol = solution.primalRowSolution();
2035#ifndef NDEBUG
2036        for (iRow = 0; iRow < numberRows; iRow++) {
2037          double value = result[iRow];
2038          assert(eq(value, sol[iRow]));
2039          assert(value < upper[iRow] + 2.0e-8);
2040          assert(value > lower[iRow] - 2.0e-8);
2041        }
2042#endif
2043        CoinFillN(result, numberRows, 0.0);
2044        solution.matrix()->times(ray, result);
2045        // there may be small differences (especially if scaled)
2046        for (iRow = 0; iRow < numberRows; iRow++) {
2047          double value = result[iRow];
2048          if (value > 1.0e-8)
2049            assert(upper[iRow] > 1.0e30);
2050          else if (value < -1.0e-8)
2051            assert(lower[iRow] < -1.0e30);
2052        }
2053        delete[] result;
2054        delete[] ray;
2055      }
2056    } else {
2057      std::cerr << "Error reading brandy from sample data. Skipping test." << std::endl;
2058    }
2059  }
2060  // test infeasible
2061  {
2062    CoinMpsIO m;
2063    std::string fn = dirSample + "brandy";
2064    if (m.readMps(fn.c_str(), "mps") == 0) {
2065      ClpSimplex solution;
2066      // do twice - without and with scaling
2067      int iPass;
2068      for (iPass = 0; iPass < 2; iPass++) {
2069        solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
2070          m.getObjCoefficients(),
2071          m.getRowLower(), m.getRowUpper());
2072        if (iPass)
2073          solution.scaling();
2074        // test infeasible and ray
2075        solution.columnUpper()[0] = 0.0;
2076#if 1 //def DUAL
2077        solution.setDualBound(100.0);
2078        solution.dual();
2079#else
2080        solution.primal();
2081#endif
2082        assert(solution.status() == 1);
2083        int numberColumns = solution.numberColumns();
2084        int numberRows = solution.numberRows();
2085        double *lower = solution.rowLower();
2086        double *upper = solution.rowUpper();
2087        double *ray = solution.infeasibilityRay();
2088        assert(ray);
2089        // construct proof of infeasibility
2090        int iRow, iColumn;
2091        double lo = 0.0, up = 0.0;
2092        int nl = 0, nu = 0;
2093        for (iRow = 0; iRow < numberRows; iRow++) {
2094          if (lower[iRow] > -1.0e20) {
2095            lo += ray[iRow] * lower[iRow];
2096          } else {
2097            if (ray[iRow] > 1.0e-8)
2098              nl++;
2099          }
2100          if (upper[iRow] < 1.0e20) {
2101            up += ray[iRow] * upper[iRow];
2102          } else {
2103            if (ray[iRow] > 1.0e-8)
2104              nu++;
2105          }
2106        }
2107        if (nl)
2108          lo = -1.0e100;
2109        if (nu)
2110          up = 1.0e100;
2111        double *result = new double[numberColumns];
2112        double lo2 = 0.0, up2 = 0.0;
2113        CoinFillN(result, numberColumns, 0.0);
2114        solution.matrix()->transposeTimes(ray, result);
2115        lower = solution.columnLower();
2116        upper = solution.columnUpper();
2117        nl = nu = 0;
2118        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2119          if (result[iColumn] > 1.0e-8) {
2120            if (lower[iColumn] > -1.0e20)
2121              lo2 += result[iColumn] * lower[iColumn];
2122            else
2123              nl++;
2124            if (upper[iColumn] < 1.0e20)
2125              up2 += result[iColumn] * upper[iColumn];
2126            else
2127              nu++;
2128          } else if (result[iColumn] < -1.0e-8) {
2129            if (lower[iColumn] > -1.0e20)
2130              up2 += result[iColumn] * lower[iColumn];
2131            else
2132              nu++;
2133            if (upper[iColumn] < 1.0e20)
2134              lo2 += result[iColumn] * upper[iColumn];
2135            else
2136              nl++;
2137          }
2138        }
2139        if (nl)
2140          lo2 = -1.0e100;
2141        if (nu)
2142          up2 = 1.0e100;
2143        // make sure inconsistency
2144        assert(lo2 > up || up2 < lo);
2145        delete[] result;
2146        delete[] ray;
2147      }
2148    } else {
2149      std::cerr << "Error reading brandy from sample data. Skipping test." << std::endl;
2150    }
2151  }
2152  // test delete and add
2153  {
2154    CoinMpsIO m;
2155    std::string fn = dirSample + "brandy";
2156    if (m.readMps(fn.c_str(), "mps") == 0) {
2157      ClpSimplex solution;
2158      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
2159        m.getObjCoefficients(),
2160        m.getRowLower(), m.getRowUpper());
2161      solution.dual();
2162      CoinRelFltEq eq(1.0e-8);
2163      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
2164
2165      int numberColumns = solution.numberColumns();
2166      int numberRows = solution.numberRows();
2167      double *saveObj = new double[numberColumns];
2168      double *saveLower = new double[numberRows + numberColumns];
2169      double *saveUpper = new double[numberRows + numberColumns];
2170      int *which = new int[numberRows + numberColumns];
2171
2172      CoinBigIndex numberElements = m.getMatrixByCol()->getNumElements();
2173      CoinBigIndex *starts = new CoinBigIndex[numberRows + numberColumns];
2174      int *index = new int[numberElements];
2175      double *element = new double[numberElements];
2176
2177      const CoinBigIndex *startM;
2178      const int *lengthM;
2179      const int *indexM;
2180      const double *elementM;
2181
2182      int n, nel;
2183
2184      // delete non basic columns
2185      n = 0;
2186      nel = 0;
2187      int iRow, iColumn;
2188      const double *lower = m.getColLower();
2189      const double *upper = m.getColUpper();
2190      const double *objective = m.getObjCoefficients();
2191      startM = m.getMatrixByCol()->getVectorStarts();
2192      lengthM = m.getMatrixByCol()->getVectorLengths();
2193      indexM = m.getMatrixByCol()->getIndices();
2194      elementM = m.getMatrixByCol()->getElements();
2195      starts[0] = 0;
2196      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2197        if (solution.getColumnStatus(iColumn) != ClpSimplex::basic) {
2198          saveObj[n] = objective[iColumn];
2199          saveLower[n] = lower[iColumn];
2200          saveUpper[n] = upper[iColumn];
2201          CoinBigIndex j;
2202          for (j = startM[iColumn]; j < startM[iColumn] + lengthM[iColumn]; j++) {
2203            index[nel] = indexM[j];
2204            element[nel++] = elementM[j];
2205          }
2206          which[n++] = iColumn;
2207          starts[n] = nel;
2208        }
2209      }
2210      solution.deleteColumns(n, which);
2211      solution.dual();
2212      // Put back
2213      solution.addColumns(n, saveLower, saveUpper, saveObj,
2214        starts, index, element);
2215      solution.dual();
2216      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
2217      // Delete all columns and add back
2218      n = 0;
2219      nel = 0;
2220      starts[0] = 0;
2221      lower = m.getColLower();
2222      upper = m.getColUpper();
2223      objective = m.getObjCoefficients();
2224      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2225        saveObj[n] = objective[iColumn];
2226        saveLower[n] = lower[iColumn];
2227        saveUpper[n] = upper[iColumn];
2228        CoinBigIndex j;
2229        for (j = startM[iColumn]; j < startM[iColumn] + lengthM[iColumn]; j++) {
2230          index[nel] = indexM[j];
2231          element[nel++] = elementM[j];
2232        }
2233        which[n++] = iColumn;
2234        starts[n] = nel;
2235      }
2236      solution.deleteColumns(n, which);
2237      solution.dual();
2238      // Put back
2239      solution.addColumns(n, saveLower, saveUpper, saveObj,
2240        starts, index, element);
2241      solution.dual();
2242      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
2243
2244      // reload with original
2245      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
2246        m.getObjCoefficients(),
2247        m.getRowLower(), m.getRowUpper());
2248      // delete half rows
2249      n = 0;
2250      nel = 0;
2251      lower = m.getRowLower();
2252      upper = m.getRowUpper();
2253      startM = m.getMatrixByRow()->getVectorStarts();
2254      lengthM = m.getMatrixByRow()->getVectorLengths();
2255      indexM = m.getMatrixByRow()->getIndices();
2256      elementM = m.getMatrixByRow()->getElements();
2257      starts[0] = 0;
2258      for (iRow = 0; iRow < numberRows; iRow++) {
2259        if ((iRow & 1) == 0) {
2260          saveLower[n] = lower[iRow];
2261          saveUpper[n] = upper[iRow];
2262          CoinBigIndex j;
2263          for (j = startM[iRow]; j < startM[iRow] + lengthM[iRow]; j++) {
2264            index[nel] = indexM[j];
2265            element[nel++] = elementM[j];
2266          }
2267          which[n++] = iRow;
2268          starts[n] = nel;
2269        }
2270      }
2271      solution.deleteRows(n, which);
2272      solution.dual();
2273      // Put back
2274      solution.addRows(n, saveLower, saveUpper,
2275        starts, index, element);
2276      solution.dual();
2277      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
2278      solution.writeMps("yy.mps");
2279      // Delete all rows
2280      n = 0;
2281      nel = 0;
2282      lower = m.getRowLower();
2283      upper = m.getRowUpper();
2284      starts[0] = 0;
2285      for (iRow = 0; iRow < numberRows; iRow++) {
2286        saveLower[n] = lower[iRow];
2287        saveUpper[n] = upper[iRow];
2288        CoinBigIndex j;
2289        for (j = startM[iRow]; j < startM[iRow] + lengthM[iRow]; j++) {
2290          index[nel] = indexM[j];
2291          element[nel++] = elementM[j];
2292        }
2293        which[n++] = iRow;
2294        starts[n] = nel;
2295      }
2296      solution.deleteRows(n, which);
2297      solution.dual();
2298      // Put back
2299      solution.addRows(n, saveLower, saveUpper,
2300        starts, index, element);
2301      solution.dual();
2302      solution.writeMps("xx.mps");
2303      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
2304      // Zero out status array to give some interest
2305      memset(solution.statusArray() + numberColumns, 0, numberRows);
2306      solution.primal(1);
2307      assert(eq(solution.objectiveValue(), 1.5185098965e+03));
2308      // Delete all columns and rows
2309      n = 0;
2310      for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2311        which[n++] = iColumn;
2312        starts[n] = nel;
2313      }
2314      solution.deleteColumns(n, which);
2315      n = 0;
2316      for (iRow = 0; iRow < numberRows; iRow++) {
2317        which[n++] = iRow;
2318        starts[n] = nel;
2319      }
2320      solution.deleteRows(n, which);
2321
2322      delete[] saveObj;
2323      delete[] saveLower;
2324      delete[] saveUpper;
2325      delete[] which;
2326      delete[] starts;
2327      delete[] index;
2328      delete[] element;
2329    } else {
2330      std::cerr << "Error reading brandy from sample data. Skipping test." << std::endl;
2331    }
2332  }
2333#if 1
2334  // Test barrier
2335  {
2336    CoinMpsIO m;
2337    std::string fn = dirSample + "exmip1";
2338    if (m.readMps(fn.c_str(), "mps") == 0) {
2339      ClpInterior solution;
2340      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
2341        m.getObjCoefficients(),
2342        m.getRowLower(), m.getRowUpper());
2343      solution.primalDual();
2344    } else {
2345      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
2346    }
2347  }
2348#endif
2349#if COIN_BIG_INDEX == 0
2350  // test network
2351#define QUADRATIC
2352  if (1) {
2353    std::string fn = dirSample + "input.130";
2354    int numberColumns;
2355    int numberRows;
2356
2357    FILE *fp = fopen(fn.c_str(), "r");
2358    if (!fp) {
2359      // Try in Data/Sample
2360      fn = "Data/Sample/input.130";
2361      fp = fopen(fn.c_str(), "r");
2362    }
2363    if (!fp) {
2364      fprintf(stderr, "Unable to open file input.130 in dirSample or Data/Sample directory\n");
2365    } else {
2366      int problem;
2367      char temp[100];
2368      // read and skip
2369      int x = fscanf(fp, "%s", temp);
2370      if (x < 0)
2371        throw("bad fscanf");
2372      assert(!strcmp(temp, "BEGIN"));
2373      x = fscanf(fp, "%*s %*s %d %d %*s %*s %d %*s", &problem, &numberRows,
2374        &numberColumns);
2375      if (x < 0)
2376        throw("bad fscanf");
2377      // scan down to SUPPLY
2378      while (fgets(temp, 100, fp)) {
2379        if (!strncmp(temp, "SUPPLY", 6))
2380          break;
2381      }
2382      if (strncmp(temp, "SUPPLY", 6)) {
2383        fprintf(stderr, "Unable to find SUPPLY\n");
2384        exit(2);
2385      }
2386      // get space for rhs
2387      double *lower = new double[numberRows];
2388      double *upper = new double[numberRows];
2389      int i;
2390      for (i = 0; i < numberRows; i++) {
2391        lower[i] = 0.0;
2392        upper[i] = 0.0;
2393      }
2394      // ***** Remember to convert to C notation
2395      while (fgets(temp, 100, fp)) {
2396        int row;
2397        int value;
2398        if (!strncmp(temp, "ARCS", 4))
2399          break;
2400        sscanf(temp, "%d %d", &row, &value);
2401        upper[row - 1] = -value;
2402        lower[row - 1] = -value;
2403      }
2404      if (strncmp(temp, "ARCS", 4)) {
2405        fprintf(stderr, "Unable to find ARCS\n");
2406        exit(2);
2407      }
2408      // number of columns may be underestimate
2409      int *head = new int[2 * numberColumns];
2410      int *tail = new int[2 * numberColumns];
2411      int *ub = new int[2 * numberColumns];
2412      int *cost = new int[2 * numberColumns];
2413      // ***** Remember to convert to C notation
2414      numberColumns = 0;
2415      while (fgets(temp, 100, fp)) {
2416        int iHead;
2417        int iTail;
2418        int iUb;
2419        int iCost;
2420        if (!strncmp(temp, "DEMAND", 6))
2421          break;
2422        sscanf(temp, "%d %d %d %d", &iHead, &iTail, &iCost, &iUb);
2423        iHead--;
2424        iTail--;
2425        head[numberColumns] = iHead;
2426        tail[numberColumns] = iTail;
2427        ub[numberColumns] = iUb;
2428        cost[numberColumns] = iCost;
2429        numberColumns++;
2430      }
2431      if (strncmp(temp, "DEMAND", 6)) {
2432        fprintf(stderr, "Unable to find DEMAND\n");
2433        exit(2);
2434      }
2435      // ***** Remember to convert to C notation
2436      while (fgets(temp, 100, fp)) {
2437        int row;
2438        int value;
2439        if (!strncmp(temp, "END", 3))
2440          break;
2441        sscanf(temp, "%d %d", &row, &value);
2442        upper[row - 1] = value;
2443        lower[row - 1] = value;
2444      }
2445      if (strncmp(temp, "END", 3)) {
2446        fprintf(stderr, "Unable to find END\n");
2447        exit(2);
2448      }
2449      printf("Problem %d has %d rows and %d columns\n", problem,
2450        numberRows, numberColumns);
2451      fclose(fp);
2452      ClpSimplex model;
2453      // now build model
2454
2455      double *objective = new double[numberColumns];
2456      double *lowerColumn = new double[numberColumns];
2457      double *upperColumn = new double[numberColumns];
2458
2459      double *element = new double[2 * numberColumns];
2460      CoinBigIndex *start = new CoinBigIndex[numberColumns + 1];
2461      int *row = new int[2 * numberColumns];
2462      start[numberColumns] = 2 * numberColumns;
2463      for (i = 0; i < numberColumns; i++) {
2464        start[i] = 2 * i;
2465        element[2 * i] = -1.0;
2466        element[2 * i + 1] = 1.0;
2467        row[2 * i] = head[i];
2468        row[2 * i + 1] = tail[i];
2469        lowerColumn[i] = 0.0;
2470        upperColumn[i] = ub[i];
2471        objective[i] = cost[i];
2472      }
2473      // Create Packed Matrix
2474      CoinPackedMatrix matrix;
2475      int *lengths = NULL;
2476      matrix.assignMatrix(true, numberRows, numberColumns,
2477        2 * numberColumns, element, row, start, lengths);
2478      // load model
2479      model.loadProblem(matrix,
2480        lowerColumn, upperColumn, objective,
2481        lower, upper);
2482      model.factorization()->maximumPivots(200 + model.numberRows() / 100);
2483      model.createStatus();
2484      double time1 = CoinCpuTime();
2485      model.dual();
2486      std::cout << "Network problem, ClpPackedMatrix took " << CoinCpuTime() - time1 << " seconds" << std::endl;
2487      ClpPlusMinusOneMatrix *plusMinus = new ClpPlusMinusOneMatrix(matrix);
2488      assert(plusMinus->getIndices()); // would be zero if not +- one
2489      //ClpPlusMinusOneMatrix *plusminus_matrix;
2490
2491      //plusminus_matrix = new ClpPlusMinusOneMatrix;
2492
2493      //plusminus_matrix->passInCopy(numberRows, numberColumns, true, plusMinus->getMutableIndices(),
2494      //                         plusMinus->startPositive(),plusMinus->startNegative());
2495      model.loadProblem(*plusMinus,
2496        lowerColumn, upperColumn, objective,
2497        lower, upper);
2498      //model.replaceMatrix( plusminus_matrix , true);
2499      delete plusMinus;
2500      //model.createStatus();
2501      //model.initialSolve();
2502      //model.writeMps("xx.mps");
2503
2504      model.factorization()->maximumPivots(200 + model.numberRows() / 100);
2505      model.createStatus();
2506      time1 = CoinCpuTime();
2507      model.dual();
2508      std::cout << "Network problem, ClpPlusMinusOneMatrix took " << CoinCpuTime() - time1 << " seconds" << std::endl;
2509      ClpNetworkMatrix network(numberColumns, head, tail);
2510      model.loadProblem(network,
2511        lowerColumn, upperColumn, objective,
2512        lower, upper);
2513
2514      model.factorization()->maximumPivots(200 + model.numberRows() / 100);
2515      model.createStatus();
2516      time1 = CoinCpuTime();
2517      model.dual();
2518      std::cout << "Network problem, ClpNetworkMatrix took " << CoinCpuTime() - time1 << " seconds" << std::endl;
2519      delete[] lower;
2520      delete[] upper;
2521      delete[] head;
2522      delete[] tail;
2523      delete[] ub;
2524      delete[] cost;
2525      delete[] objective;
2526      delete[] lowerColumn;
2527      delete[] upperColumn;
2528    }
2529  }
2530#ifdef QUADRATIC
2531  // Test quadratic to solve linear
2532  if (1) {
2533    CoinMpsIO m;
2534    std::string fn = dirSample + "exmip1";
2535    if (m.readMps(fn.c_str(), "mps") == 0) {
2536      ClpSimplex solution;
2537      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
2538        m.getObjCoefficients(),
2539        m.getRowLower(), m.getRowUpper());
2540      //solution.dual();
2541      // get quadratic part
2542      int numberColumns = solution.numberColumns();
2543      int *start = new int[numberColumns + 1];
2544      int *column = new int[numberColumns];
2545      double *element = new double[numberColumns];
2546      int i;
2547      start[0] = 0;
2548      int n = 0;
2549      int kk = numberColumns - 1;
2550      int kk2 = numberColumns - 1;
2551      for (i = 0; i < numberColumns; i++) {
2552        if (i >= kk) {
2553          column[n] = i;
2554          if (i >= kk2)
2555            element[n] = 1.0e-1;
2556          else
2557            element[n] = 0.0;
2558          n++;
2559        }
2560        start[i + 1] = n;
2561      }
2562      // Load up objective
2563      solution.loadQuadraticObjective(numberColumns, start, column, element);
2564      delete[] start;
2565      delete[] column;
2566      delete[] element;
2567      //solution.quadraticSLP(50,1.0e-4);
2568      CoinRelFltEq eq(1.0e-4);
2569      //assert(eq(objValue,820.0));
2570      //solution.setLogLevel(63);
2571      solution.primal();
2572      printSol(solution);
2573      //assert(eq(objValue,3.2368421));
2574      //exit(77);
2575    } else {
2576      std::cerr << "Error reading exmip1 from sample data. Skipping test." << std::endl;
2577    }
2578  }
2579  // Test quadratic
2580  if (1) {
2581    CoinMpsIO m;
2582    std::string fn = dirSample + "share2qp";
2583    //fn = "share2qpb";
2584    if (m.readMps(fn.c_str(), "mps") == 0) {
2585      ClpSimplex model;
2586      model.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
2587        m.getObjCoefficients(),
2588        m.getRowLower(), m.getRowUpper());
2589      model.dual();
2590      // get quadratic part
2591      int *start = NULL;
2592      int *column = NULL;
2593      double *element = NULL;
2594      m.readQuadraticMps(NULL, start, column, element, 2);
2595      int column2[200];
2596      double element2[200];
2597      int start2[80];
2598      int j;
2599      start2[0] = 0;
2600      int nel = 0;
2601      bool good = false;
2602      for (j = 0; j < 79; j++) {
2603        if (start[j] == start[j + 1]) {
2604          column2[nel] = j;
2605          element2[nel] = 0.0;
2606          nel++;
2607        } else {
2608          int i;
2609          for (i = start[j]; i < start[j + 1]; i++) {
2610            column2[nel] = column[i];
2611            element2[nel++] = element[i];
2612          }
2613        }
2614        start2[j + 1] = nel;
2615      }
2616      // Load up objective
2617      if (good)
2618        model.loadQuadraticObjective(model.numberColumns(), start2, column2, element2);
2619      else
2620        model.loadQuadraticObjective(model.numberColumns(), start, column, element);
2621      delete[] start;
2622      delete[] column;
2623      delete[] element;
2624      int numberColumns = model.numberColumns();
2625      model.scaling(0);
2626#if 0
2627               model.nonlinearSLP(50, 1.0e-4);
2628#else
2629      // Get feasible
2630      ClpObjective *saveObjective = model.objectiveAsObject()->clone();
2631      ClpLinearObjective zeroObjective(NULL, numberColumns);
2632      model.setObjective(&zeroObjective);
2633      model.dual();
2634      model.setObjective(saveObjective);
2635      delete saveObjective;
2636#endif
2637      //model.setLogLevel(63);
2638      //exit(77);
2639      model.setFactorizationFrequency(10);
2640      model.primal();
2641      printSol(model);
2642#ifndef NDEBUG
2643      double objValue = model.getObjValue();
2644#endif
2645      CoinRelFltEq eq(1.0e-4);
2646      assert(eq(objValue, -400.92));
2647      // and again for barrier
2648      model.barrier(false);
2649      //printSol(model);
2650      model.allSlackBasis();
2651      model.primal();
2652      //printSol(model);
2653    } else {
2654      std::cerr << "Error reading share2qp from sample data. Skipping test." << std::endl;
2655    }
2656  }
2657  if (0) {
2658    CoinMpsIO m;
2659    std::string fn = "./beale";
2660    //fn = "./jensen";
2661    if (m.readMps(fn.c_str(), "mps") == 0) {
2662      ClpSimplex solution;
2663      solution.loadProblem(*m.getMatrixByCol(), m.getColLower(), m.getColUpper(),
2664        m.getObjCoefficients(),
2665        m.getRowLower(), m.getRowUpper());
2666      solution.setDblParam(ClpObjOffset, m.objectiveOffset());
2667      solution.dual();
2668      // get quadratic part
2669      int *start = NULL;
2670      int *column = NULL;
2671      double *element = NULL;
2672      m.readQuadraticMps(NULL, start, column, element, 2);
2673      // Load up objective
2674      solution.loadQuadraticObjective(solution.numberColumns(), start, column, element);
2675      delete[] start;
2676      delete[] column;
2677      delete[] element;
2678      solution.primal(1);
2679      solution.nonlinearSLP(50, 1.0e-4);
2680      double objValue = solution.getObjValue();
2681      CoinRelFltEq eq(1.0e-4);
2682      assert(eq(objValue, 0.5));
2683      solution.primal();
2684      objValue = solution.getObjValue();
2685      assert(eq(objValue, 0.5));
2686    } else {
2687      std::cerr << "Error reading beale.mps. Skipping test." << std::endl;
2688    }
2689  }
2690#endif
2691  // Test CoinStructuredModel
2692  {
2693
2694    // Sub block
2695    CoinModel sub;
2696    {
2697      // matrix data
2698      //deliberate hiccup of 2 between 0 and 1
2699      CoinBigIndex start[5] = { 0, 4, 7, 8, 9 };
2700      int length[5] = { 2, 3, 1, 1, 1 };
2701      int rows[11] = { 0, 2, -1, -1, 0, 1, 2, 0, 1, 2 };
2702      double elements[11] = { 7.0, 2.0, 1.0e10, 1.0e10, -2.0, 1.0, -2.0, 1, 1, 1 };
2703      CoinPackedMatrix matrix(true, 3, 5, 8, elements, rows, start, length);
2704      // by row
2705      matrix.reverseOrdering();
2706      const double *element = matrix.getElements();
2707      const int *column = matrix.getIndices();
2708      const CoinBigIndex *rowStart = matrix.getVectorStarts();
2709      const int *rowLength = matrix.getVectorLengths();
2710
2711      // rim data
2712      //double objective[5]={-4.0,1.0,0.0,0.0,0.0};
2713      double rowLower[3] = { 14.0, 3.0, 3.0 };
2714      double rowUpper[3] = { 14.0, 3.0, 3.0 };
2715      //double columnLower[5]={0.0,0.0,0.0,0.0,0.0};
2716      //double columnUpper[5]={100.0,100.0,100.0,100.0,100.0};
2717
2718      for (int i = 0; i < 3; i++) {
2719        sub.addRow(rowLength[i], column + rowStart[i],
2720          element + rowStart[i], rowLower[i], rowUpper[i]);
2721      }
2722      //for (int i=0;i<5;i++) {
2723      //sub.setColumnBounds(i,columnLower[i],columnUpper[i]);
2724      //sub.setColumnObjective(i,objective[i]);
2725      //}
2726      sub.convertMatrix();
2727    }
2728    // Top
2729    CoinModel top;
2730    {
2731      // matrix data
2732      CoinBigIndex start[5] = { 0, 2, 4, 6, 8 };
2733      int length[5] = { 2, 2, 2, 2, 2 };
2734      int rows[10] = { 0, 1, 0, 1, 0, 1, 0, 1, 0, 1 };
2735      double elements[10] = { 1, 1, 1, 1, 1, 1, 1, 1, 1, 1 };
2736      CoinPackedMatrix matrix(true, 2, 5, 8, elements, rows, start, length);
2737      // by row
2738      matrix.reverseOrdering();
2739      const double *element = matrix.getElements();
2740      const int *column = matrix.getIndices();
2741      const CoinBigIndex *rowStart = matrix.getVectorStarts();
2742      const int *rowLength = matrix.getVectorLengths();
2743
2744      // rim data
2745      double objective[5] = { -4.0, 1.0, 0.0, 0.0, 0.0 };
2746      //double rowLower[3]={14.0,3.0,3.0};
2747      //double rowUpper[3]={14.0,3.0,3.0};
2748      double columnLower[5] = { 0.0, 0.0, 0.0, 0.0, 0.0 };
2749      double columnUpper[5] = { 100.0, 100.0, 100.0, 100.0, 100.0 };
2750
2751      for (int i = 0; i < 2; i++) {
2752        top.addRow(rowLength[i], column + rowStart[i],
2753          element + rowStart[i],
2754          -COIN_DBL_MAX, COIN_DBL_MAX);
2755      }
2756      for (int i = 0; i < 5; i++) {
2757        top.setColumnBounds(i, columnLower[i], columnUpper[i]);
2758        top.setColumnObjective(i, objective[i]);
2759      }
2760      top.convertMatrix();
2761    }
2762    // Create a structured model
2763    CoinStructuredModel structured;
2764    int numberBlocks = 5;
2765    for (int i = 0; i < numberBlocks; i++) {
2766      std::string topName = "row_master";
2767      std::string blockName = "block_";
2768      char bName = static_cast< char >('a' + static_cast< char >(i));
2769      blockName.append(1, bName);
2770      structured.addBlock(topName, blockName, top);
2771      structured.addBlock(blockName, blockName, sub);
2772    }
2773    // Set rhs on first block
2774    CoinModel *first = structured.coinBlock(0);
2775    for (int i = 0; i < 2; i++) {
2776      first->setRowLower(i, 0.0);
2777      first->setRowUpper(i, 100.0);
2778    }
2779    // Refresh whats set
2780    structured.refresh(0);
2781    // Could perturb stuff, but for first go don't bother
2782    ClpSimplex fullModel;
2783    // There is no original stuff set - think
2784    fullModel.loadProblem(structured, false);
2785    fullModel.dual();
2786    fullModel.dropNames();
2787    fullModel.writeMps("test.mps");
2788    // Make up very simple nested model - not realistic
2789    // Create a structured model
2790    CoinStructuredModel structured2;
2791    numberBlocks = 3;
2792    for (int i = 0; i < numberBlocks; i++) {
2793      std::string blockName = "block_";
2794      char bName = static_cast< char >('a' + static_cast< char >(i));
2795      blockName.append(1, bName);
2796      structured2.addBlock(blockName, blockName, structured);
2797    }
2798    fullModel.loadProblem(structured2, false);
2799    fullModel.dual();
2800    fullModel.dropNames();
2801    fullModel.writeMps("test2.mps");
2802  }
2803#endif
2804}
2805
2806/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
2807*/
Note: See TracBrowser for help on using the repository browser.