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

Last change on this file since 1559 was 1559, checked in by stefan, 10 years ago

merge split branch into trunk

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