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

Last change on this file since 1726 was 1726, checked in by stefan, 8 years ago

fix compiler warnings, including one that pointed to a bug

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