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

Last change on this file since 1693 was 1693, checked in by stefan, 9 years ago

make WSMP recoginition and use similar to Ipopt's WSMP interface

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