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

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

minor changes to implement Aboca

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