Context Navigation

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

Last change on this file since 1693 was 1693, checked in by stefan, 9 years ago
make WSMP recoginition and use similar to Ipopt's WSMP interface
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 107.9 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: