Context Navigation

source: trunk/Clp/src/ClpMain.cpp @ 1344

Last change on this file since 1344 was 1344, checked in by forrest, 12 years ago
changes to simplex and lots of stuff and start Mumps cholesky
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 91.4 KB

Line
1	// Copyright (C) 2002, International Business Machines
2	// Corporation and others. All Rights Reserved.
3
4	#include "CoinPragma.hpp"
5
6	#include <cassert>
7	#include <cstdio>
8	#include <cstdlib>
9	#include <cmath>
10	#include <cfloat>
11	#include <string>
12	#include <iostream>
13
14	int boundary_sort=1000;
15	int boundary_sort2=1000;
16	int boundary_sort3=10000;
17
18	#include "CoinPragma.hpp"
19	#include "CoinHelperFunctions.hpp"
20	#include "CoinSort.hpp"
21	// History since 1.0 at end
22	#define CLPVERSION "1.09.00"
23
24	#include "CoinMpsIO.hpp"
25	#include "CoinFileIO.hpp"
26
27	#include "ClpFactorization.hpp"
28	#include "CoinTime.hpp"
29	#include "ClpSimplex.hpp"
30	#include "ClpSimplexOther.hpp"
31	#include "ClpSolve.hpp"
32	#include "ClpPackedMatrix.hpp"
33	#include "ClpPlusMinusOneMatrix.hpp"
34	#include "ClpNetworkMatrix.hpp"
35	#include "ClpDualRowSteepest.hpp"
36	#include "ClpDualRowDantzig.hpp"
37	#include "ClpLinearObjective.hpp"
38	#include "ClpPrimalColumnSteepest.hpp"
39	#include "ClpPrimalColumnDantzig.hpp"
40	#include "ClpPresolve.hpp"
41	#include "CbcOrClpParam.hpp"
42	#include "CoinSignal.hpp"
43	#ifdef DMALLOC
44	#include "dmalloc.h"
45	#endif
46	#ifdef WSSMP_BARRIER
47	#define FOREIGN_BARRIER
48	#endif
49	#ifdef UFL_BARRIER
50	#define FOREIGN_BARRIER
51	#endif
52	#ifdef TAUCS_BARRIER
53	#define FOREIGN_BARRIER
54	#endif
55	#ifdef MUMPS_BARRIER
56	#define FOREIGN_BARRIER
57	#endif
58
59	static double totalTime=0.0;
60	static bool maskMatches(const int * starts, char ** masks,
61	std::string & check);
62	static ClpSimplex * currentModel = NULL;
63
64	extern "C" {
65	static void
66	#if defined(_MSC_VER)
67	__cdecl
68	#endif // _MSC_VER
69	signal_handler(int whichSignal)
70	{
71	if (currentModel!=NULL)
72	currentModel->setMaximumIterations(0); // stop at next iterations
73	return;
74	}
75	}
76
77	//#############################################################################
78
79	#ifdef NDEBUG
80	#undef NDEBUG
81	#endif
82
83	int mainTest (int argc, const char *argv[],int algorithm,
84	ClpSimplex empty, bool doPresolve,int switchOff,bool doVector);
85	static void statistics(ClpSimplex * originalModel, ClpSimplex * model);
86	static void generateCode(const char * fileName,int type);
87	// Returns next valid field
88	int CbcOrClpRead_mode=1;
89	FILE * CbcOrClpReadCommand=stdin;
90
91	int
92	#if defined(_MSC_VER)
93	__cdecl
94	#endif // _MSC_VER
95	main (int argc, const char *argv[])
96	{
97	// next {} is just to make sure all memory should be freed - for debug
98	{
99	double time1 = CoinCpuTime(),time2;
100	// Set up all non-standard stuff
101	//int numberModels=1;
102	ClpSimplex * models = new ClpSimplex[1];
103
104	// default action on import
105	int allowImportErrors=0;
106	int keepImportNames=1;
107	int doIdiot=-1;
108	int outputFormat=2;
109	int slpValue=-1;
110	int cppValue=-1;
111	int printOptions=0;
112	int printMode=0;
113	int presolveOptions=0;
114	int doCrash=0;
115	int doVector=0;
116	int doSprint=-1;
117	// set reasonable defaults
118	int preSolve=5;
119	bool preSolveFile=false;
120	models->setPerturbation(50);
121	models->messageHandler()->setPrefix(false);
122	const char dirsep = CoinFindDirSeparator();
123	std::string directory;
124	std::string dirSample;
125	std::string dirNetlib;
126	std::string dirMiplib;
127	if (dirsep == '/') {
128	directory = "./";
129	dirSample = "../../Data/Sample/";
130	dirNetlib = "../../Data/Netlib/";
131	dirMiplib = "../../Data/miplib3/";
132	} else {
133	directory = ".\\";
134	dirSample = "..\\..\\Data\\Sample\\";
135	dirNetlib = "..\\..\\Data\\Netlib\\";
136	dirMiplib = "..\\..\\Data\\miplib3\\";
137	}
138	std::string defaultDirectory = directory;
139	std::string importFile ="";
140	std::string exportFile ="default.mps";
141	std::string importBasisFile ="";
142	int basisHasValues=0;
143	int substitution=3;
144	int dualize=3; // dualize if looks promising
145	std::string exportBasisFile ="default.bas";
146	std::string saveFile ="default.prob";
147	std::string restoreFile ="default.prob";
148	std::string solutionFile ="stdout";
149	std::string solutionSaveFile ="solution.file";
150	std::string printMask="";
151	CbcOrClpParam parameters[CBCMAXPARAMETERS];
152	int numberParameters ;
153	establishParams(numberParameters,parameters) ;
154	parameters[whichParam(BASISIN,numberParameters,parameters)].setStringValue(importBasisFile);
155	parameters[whichParam(BASISOUT,numberParameters,parameters)].setStringValue(exportBasisFile);
156	parameters[whichParam(PRINTMASK,numberParameters,parameters)].setStringValue(printMask);
157	parameters[whichParam(DIRECTORY,numberParameters,parameters)].setStringValue(directory);
158	parameters[whichParam(DIRSAMPLE,numberParameters,parameters)].setStringValue(dirSample);
159	parameters[whichParam(DIRNETLIB,numberParameters,parameters)].setStringValue(dirNetlib);
160	parameters[whichParam(DIRMIPLIB,numberParameters,parameters)].setStringValue(dirMiplib);
161	parameters[whichParam(DUALBOUND,numberParameters,parameters)].setDoubleValue(models->dualBound());
162	parameters[whichParam(DUALTOLERANCE,numberParameters,parameters)].setDoubleValue(models->dualTolerance());
163	parameters[whichParam(EXPORT,numberParameters,parameters)].setStringValue(exportFile);
164	parameters[whichParam(IDIOT,numberParameters,parameters)].setIntValue(doIdiot);
165	parameters[whichParam(IMPORT,numberParameters,parameters)].setStringValue(importFile);
166	parameters[whichParam(SOLVERLOGLEVEL,numberParameters,parameters)].setIntValue(models->logLevel());
167	parameters[whichParam(MAXFACTOR,numberParameters,parameters)].setIntValue(models->factorizationFrequency());
168	parameters[whichParam(MAXITERATION,numberParameters,parameters)].setIntValue(models->maximumIterations());
169	parameters[whichParam(OUTPUTFORMAT,numberParameters,parameters)].setIntValue(outputFormat);
170	parameters[whichParam(PRESOLVEPASS,numberParameters,parameters)].setIntValue(preSolve);
171	parameters[whichParam(PERTVALUE,numberParameters,parameters)].setIntValue(models->perturbation());
172	parameters[whichParam(PRIMALTOLERANCE,numberParameters,parameters)].setDoubleValue(models->primalTolerance());
173	parameters[whichParam(PRIMALWEIGHT,numberParameters,parameters)].setDoubleValue(models->infeasibilityCost());
174	parameters[whichParam(RESTORE,numberParameters,parameters)].setStringValue(restoreFile);
175	parameters[whichParam(SAVE,numberParameters,parameters)].setStringValue(saveFile);
176	parameters[whichParam(TIMELIMIT,numberParameters,parameters)].setDoubleValue(models->maximumSeconds());
177	parameters[whichParam(SOLUTION,numberParameters,parameters)].setStringValue(solutionFile);
178	parameters[whichParam(SAVESOL,numberParameters,parameters)].setStringValue(solutionSaveFile);
179	parameters[whichParam(SPRINT,numberParameters,parameters)].setIntValue(doSprint);
180	parameters[whichParam(SUBSTITUTION,numberParameters,parameters)].setIntValue(substitution);
181	parameters[whichParam(DUALIZE,numberParameters,parameters)].setIntValue(dualize);
182	parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].setDoubleValue(1.0e-8);
183	int verbose=0;
184
185	// total number of commands read
186	int numberGoodCommands=0;
187	bool * goodModels = new bool[1];
188
189	// Hidden stuff for barrier
190	int choleskyType = 0;
191	int gamma=0;
192	int scaleBarrier=0;
193	int doKKT=0;
194	int crossover=2; // do crossover unless quadratic
195
196	int iModel=0;
197	goodModels[0]=false;
198	//models[0].scaling(1);
199	//models[0].setDualBound(1.0e6);
200	//models[0].setDualTolerance(1.0e-7);
201	//ClpDualRowSteepest steep;
202	//models[0].setDualRowPivotAlgorithm(steep);
203	//models[0].setPrimalTolerance(1.0e-7);
204	//ClpPrimalColumnSteepest steepP;
205	//models[0].setPrimalColumnPivotAlgorithm(steepP);
206	std::string field;
207	std::cout<<"Coin LP version "<<CLPVERSION
208	<<", build "<<__DATE__<<std::endl;
209	// Print command line
210	if (argc>1) {
211	printf("command line - ");
212	for (int i=0;i<argc;i++)
213	printf("%s ",argv[i]);
214	printf("\n");
215	}
216
217	while (1) {
218	// next command
219	field=CoinReadGetCommand(argc,argv);
220
221	// exit if null or similar
222	if (!field.length()) {
223	if (numberGoodCommands==1&&goodModels[0]) {
224	// we just had file name - do dual or primal
225	field="either";
226	} else if (!numberGoodCommands) {
227	// let's give the sucker a hint
228	std::cout
229	<<"Clp takes input from arguments ( - switches to stdin)"
230	<<std::endl
231	<<"Enter ? for list of commands or help"<<std::endl;
232	field="-";
233	} else {
234	break;
235	}
236	}
237
238	// see if ? at end
239	int numberQuery=0;
240	if (field!="?"&&field!="???") {
241	int length = field.length();
242	int i;
243	for (i=length-1;i>0;i--) {
244	if (field[i]=='?')
245	numberQuery++;
246	else
247	break;
248	}
249	field=field.substr(0,length-numberQuery);
250	}
251	// find out if valid command
252	int iParam;
253	int numberMatches=0;
254	int firstMatch=-1;
255	for ( iParam=0; iParam<numberParameters; iParam++ ) {
256	int match = parameters[iParam].matches(field);
257	if (match==1) {
258	numberMatches = 1;
259	firstMatch=iParam;
260	break;
261	} else {
262	if (match&&firstMatch<0)
263	firstMatch=iParam;
264	numberMatches += match>>1;
265	}
266	}
267	if (iParam<numberParameters&&!numberQuery) {
268	// found
269	CbcOrClpParam found = parameters[iParam];
270	CbcOrClpParameterType type = found.type();
271	int valid;
272	numberGoodCommands++;
273	if (type==GENERALQUERY) {
274	std::cout<<"In argument list keywords have leading - "
275	", -stdin or just - switches to stdin"<<std::endl;
276	std::cout<<"One command per line (and no -)"<<std::endl;
277	std::cout<<"abcd? gives list of possibilities, if only one + explanation"<<std::endl;
278	std::cout<<"abcd?? adds explanation, if only one fuller help"<<std::endl;
279	std::cout<<"abcd without value (where expected) gives current value"<<std::endl;
280	std::cout<<"abcd value sets value"<<std::endl;
281	std::cout<<"Commands are:"<<std::endl;
282	int maxAcross=5;
283	bool evenHidden=false;
284	if ((verbose&8)!=0) {
285	// even hidden
286	evenHidden = true;
287	verbose &= ~8;
288	}
289	if (verbose)
290	maxAcross=1;
291	int limits[]={1,101,201,301,401};
292	std::vector<std::string> types;
293	types.push_back("Double parameters:");
294	types.push_back("Int parameters:");
295	types.push_back("Keyword parameters:");
296	types.push_back("Actions or string parameters:");
297	int iType;
298	for (iType=0;iType<4;iType++) {
299	int across=0;
300	if ((verbose%4)!=0)
301	std::cout<<std::endl;
302	std::cout<<types[iType]<<std::endl;
303	if ((verbose&2)!=0)
304	std::cout<<std::endl;
305	for ( iParam=0; iParam<numberParameters; iParam++ ) {
306	int type = parameters[iParam].type();
307	if ((parameters[iParam].displayThis()\|\|evenHidden)&&
308	type>=limits[iType]
309	&&type<limits[iType+1]) {
310	if (!across) {
311	if ((verbose&2)==0)
312	std::cout<<" ";
313	else
314	std::cout<<"Command ";
315	}
316	std::cout<<parameters[iParam].matchName()<<" ";
317	across++;
318	if (across==maxAcross) {
319	across=0;
320	if (verbose) {
321	// put out description as well
322	if ((verbose&1)!=0)
323	std::cout<<parameters[iParam].shortHelp();
324	std::cout<<std::endl;
325	if ((verbose&2)!=0) {
326	std::cout<<"---- description"<<std::endl;
327	parameters[iParam].printLongHelp();
328	std::cout<<"----"<<std::endl<<std::endl;
329	}
330	} else {
331	std::cout<<std::endl;
332	}
333	}
334	}
335	}
336	if (across)
337	std::cout<<std::endl;
338	}
339	} else if (type==FULLGENERALQUERY) {
340	std::cout<<"Full list of commands is:"<<std::endl;
341	int maxAcross=5;
342	int limits[]={1,101,201,301,401};
343	std::vector<std::string> types;
344	types.push_back("Double parameters:");
345	types.push_back("Int parameters:");
346	types.push_back("Keyword parameters and others:");
347	types.push_back("Actions:");
348	int iType;
349	for (iType=0;iType<4;iType++) {
350	int across=0;
351	std::cout<<types[iType]<<std::endl;
352	for ( iParam=0; iParam<numberParameters; iParam++ ) {
353	int type = parameters[iParam].type();
354	if (type>=limits[iType]
355	&&type<limits[iType+1]) {
356	if (!across)
357	std::cout<<" ";
358	std::cout<<parameters[iParam].matchName()<<" ";
359	across++;
360	if (across==maxAcross) {
361	std::cout<<std::endl;
362	across=0;
363	}
364	}
365	}
366	if (across)
367	std::cout<<std::endl;
368	}
369	} else if (type<101) {
370	// get next field as double
371	double value = CoinReadGetDoubleField(argc,argv,&valid);
372	if (!valid) {
373	parameters[iParam].setDoubleParameter(models+iModel,value);
374	} else if (valid==1) {
375	std::cout<<" is illegal for double parameter "<<parameters[iParam].name()<<" value remains "<<
376	parameters[iParam].doubleValue()<<std::endl;
377	} else {
378	std::cout<<parameters[iParam].name()<<" has value "<<
379	parameters[iParam].doubleValue()<<std::endl;
380	}
381	} else if (type<201) {
382	// get next field as int
383	int value = CoinReadGetIntField(argc,argv,&valid);
384	if (!valid) {
385	if (parameters[iParam].type()==PRESOLVEPASS)
386	preSolve = value;
387	else if (parameters[iParam].type()==IDIOT)
388	doIdiot = value;
389	else if (parameters[iParam].type()==SPRINT)
390	doSprint = value;
391	else if (parameters[iParam].type()==OUTPUTFORMAT)
392	outputFormat = value;
393	else if (parameters[iParam].type()==SLPVALUE)
394	slpValue = value;
395	else if (parameters[iParam].type()==CPP)
396	cppValue = value;
397	else if (parameters[iParam].type()==PRESOLVEOPTIONS)
398	presolveOptions = value;
399	else if (parameters[iParam].type()==PRINTOPTIONS)
400	printOptions = value;
401	else if (parameters[iParam].type()==SUBSTITUTION)
402	substitution = value;
403	else if (parameters[iParam].type()==DUALIZE)
404	dualize = value;
405	else if (parameters[iParam].type()==VERBOSE)
406	verbose = value;
407	parameters[iParam].setIntParameter(models+iModel,value);
408	} else if (valid==1) {
409	std::cout<<" is illegal for integer parameter "<<parameters[iParam].name()<<" value remains "<<
410	parameters[iParam].intValue()<<std::endl;
411	} else {
412	std::cout<<parameters[iParam].name()<<" has value "<<
413	parameters[iParam].intValue()<<std::endl;
414	}
415	} else if (type<301) {
416	// one of several strings
417	std::string value = CoinReadGetString(argc,argv);
418	int action = parameters[iParam].parameterOption(value);
419	if (action<0) {
420	if (value!="EOL") {
421	// no match
422	parameters[iParam].printOptions();
423	} else {
424	// print current value
425	std::cout<<parameters[iParam].name()<<" has value "<<
426	parameters[iParam].currentOption()<<std::endl;
427	}
428	} else {
429	parameters[iParam].setCurrentOption(action);
430	// for now hard wired
431	switch (type) {
432	case DIRECTION:
433	if (action==0)
434	models[iModel].setOptimizationDirection(1);
435	else if (action==1)
436	models[iModel].setOptimizationDirection(-1);
437	else
438	models[iModel].setOptimizationDirection(0);
439	break;
440	case DUALPIVOT:
441	if (action==0) {
442	ClpDualRowSteepest steep(3);
443	models[iModel].setDualRowPivotAlgorithm(steep);
444	} else if (action==1) {
445	//ClpDualRowDantzig dantzig;
446	ClpDualRowSteepest dantzig(5);
447	models[iModel].setDualRowPivotAlgorithm(dantzig);
448	} else if (action==2) {
449	// partial steep
450	ClpDualRowSteepest steep(2);
451	models[iModel].setDualRowPivotAlgorithm(steep);
452	} else {
453	ClpDualRowSteepest steep;
454	models[iModel].setDualRowPivotAlgorithm(steep);
455	}
456	break;
457	case PRIMALPIVOT:
458	if (action==0) {
459	ClpPrimalColumnSteepest steep(3);
460	models[iModel].setPrimalColumnPivotAlgorithm(steep);
461	} else if (action==1) {
462	ClpPrimalColumnSteepest steep(0);
463	models[iModel].setPrimalColumnPivotAlgorithm(steep);
464	} else if (action==2) {
465	ClpPrimalColumnDantzig dantzig;
466	models[iModel].setPrimalColumnPivotAlgorithm(dantzig);
467	} else if (action==3) {
468	ClpPrimalColumnSteepest steep(4);
469	models[iModel].setPrimalColumnPivotAlgorithm(steep);
470	} else if (action==4) {
471	ClpPrimalColumnSteepest steep(1);
472	models[iModel].setPrimalColumnPivotAlgorithm(steep);
473	} else if (action==5) {
474	ClpPrimalColumnSteepest steep(2);
475	models[iModel].setPrimalColumnPivotAlgorithm(steep);
476	} else if (action==6) {
477	ClpPrimalColumnSteepest steep(10);
478	models[iModel].setPrimalColumnPivotAlgorithm(steep);
479	}
480	break;
481	case SCALING:
482	models[iModel].scaling(action);
483	break;
484	case AUTOSCALE:
485	models[iModel].setAutomaticScaling(action!=0);
486	break;
487	case SPARSEFACTOR:
488	models[iModel].setSparseFactorization((1-action)!=0);
489	break;
490	case BIASLU:
491	models[iModel].factorization()->setBiasLU(action);
492	break;
493	case PERTURBATION:
494	if (action==0)
495	models[iModel].setPerturbation(50);
496	else
497	models[iModel].setPerturbation(100);
498	break;
499	case ERRORSALLOWED:
500	allowImportErrors = action;
501	break;
502	case INTPRINT:
503	printMode=action;
504	break;
505	case KEEPNAMES:
506	keepImportNames = 1-action;
507	break;
508	case PRESOLVE:
509	if (action==0)
510	preSolve = 5;
511	else if (action==1)
512	preSolve=0;
513	else if (action==2)
514	preSolve=10;
515	else
516	preSolveFile=true;
517	break;
518	case PFI:
519	models[iModel].factorization()->setForrestTomlin(action==0);
520	break;
521	case CRASH:
522	doCrash=action;
523	break;
524	case VECTOR:
525	doVector=action;
526	break;
527	case MESSAGES:
528	models[iModel].messageHandler()->setPrefix(action!=0);
529	break;
530	case CHOLESKY:
531	choleskyType = action;
532	break;
533	case GAMMA:
534	gamma=action;
535	break;
536	case BARRIERSCALE:
537	scaleBarrier=action;
538	break;
539	case KKT:
540	doKKT=action;
541	break;
542	case CROSSOVER:
543	crossover=action;
544	break;
545	default:
546	abort();
547	}
548	}
549	} else {
550	// action
551	if (type==EXIT)
552	break; // stop all
553	switch (type) {
554	case DUALSIMPLEX:
555	case PRIMALSIMPLEX:
556	case EITHERSIMPLEX:
557	case BARRIER:
558	// synonym for dual
559	case BAB:
560	if (goodModels[iModel]) {
561	double objScale =
562	parameters[whichParam(OBJSCALE2,numberParameters,parameters)].doubleValue();
563	if (objScale!=1.0) {
564	int iColumn;
565	int numberColumns=models[iModel].numberColumns();
566	double * dualColumnSolution =
567	models[iModel].dualColumnSolution();
568	ClpObjective * obj = models[iModel].objectiveAsObject();
569	assert(dynamic_cast<ClpLinearObjective *> (obj));
570	double offset;
571	double * objective = obj->gradient(NULL,NULL,offset,true);
572	for (iColumn=0;iColumn<numberColumns;iColumn++) {
573	dualColumnSolution[iColumn] *= objScale;
574	objective[iColumn] *= objScale;;
575	}
576	int iRow;
577	int numberRows=models[iModel].numberRows();
578	double * dualRowSolution =
579	models[iModel].dualRowSolution();
580	for (iRow=0;iRow<numberRows;iRow++)
581	dualRowSolution[iRow] *= objScale;
582	models[iModel].setObjectiveOffset(objScale*models[iModel].objectiveOffset());
583	}
584	ClpSolve::SolveType method;
585	ClpSolve::PresolveType presolveType;
586	ClpSimplex * model2 = models+iModel;
587	if (dualize) {
588	bool tryIt=true;
589	double fractionColumn=1.0;
590	double fractionRow=1.0;
591	if (dualize==3) {
592	dualize=1;
593	int numberColumns=model2->numberColumns();
594	int numberRows=model2->numberRows();
595	if (numberRows<50000\|\|5*numberColumns>numberRows) {
596	tryIt=false;
597	} else {
598	fractionColumn=0.1;
599	fractionRow=0.1;
600	}
601	}
602	if (tryIt) {
603	model2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel(fractionRow,fractionColumn);
604	if (model2) {
605	printf("Dual of model has %d rows and %d columns\n",
606	model2->numberRows(),model2->numberColumns());
607	model2->setOptimizationDirection(1.0);
608	} else {
609	model2 = models+iModel;
610	dualize=0;
611	}
612	} else {
613	dualize=0;
614	}
615	}
616	ClpSolve solveOptions;
617	solveOptions.setPresolveActions(presolveOptions);
618	solveOptions.setSubstitution(substitution);
619	if (preSolve!=5&&preSolve) {
620	presolveType=ClpSolve::presolveNumber;
621	if (preSolve<0) {
622	preSolve = - preSolve;
623	if (preSolve<=100) {
624	presolveType=ClpSolve::presolveNumber;
625	printf("Doing %d presolve passes - picking up non-costed slacks\n",
626	preSolve);
627	solveOptions.setDoSingletonColumn(true);
628	} else {
629	preSolve -=100;
630	presolveType=ClpSolve::presolveNumberCost;
631	printf("Doing %d presolve passes - picking up costed slacks\n",
632	preSolve);
633	}
634	}
635	} else if (preSolve) {
636	presolveType=ClpSolve::presolveOn;
637	} else {
638	presolveType=ClpSolve::presolveOff;
639	}
640	solveOptions.setPresolveType(presolveType,preSolve);
641	if (type==DUALSIMPLEX\|\|type==BAB) {
642	method=ClpSolve::useDual;
643	} else if (type==PRIMALSIMPLEX) {
644	method=ClpSolve::usePrimalorSprint;
645	} else if (type==EITHERSIMPLEX) {
646	method=ClpSolve::automatic;
647	} else {
648	method = ClpSolve::useBarrier;
649	if (crossover==1) {
650	method=ClpSolve::useBarrierNoCross;
651	} else if (crossover==2) {
652	ClpObjective * obj = models[iModel].objectiveAsObject();
653	if (obj->type()>1) {
654	method=ClpSolve::useBarrierNoCross;
655	presolveType=ClpSolve::presolveOff;
656	solveOptions.setPresolveType(presolveType,preSolve);
657	}
658	}
659	}
660	solveOptions.setSolveType(method);
661	if(preSolveFile)
662	presolveOptions \|= 0x40000000;
663	solveOptions.setSpecialOption(4,presolveOptions);
664	solveOptions.setSpecialOption(5,printOptions&1);
665	if (doVector) {
666	ClpMatrixBase * matrix = models[iModel].clpMatrix();
667	if (dynamic_cast< ClpPackedMatrix*>(matrix)) {
668	ClpPackedMatrix * clpMatrix = dynamic_cast< ClpPackedMatrix*>(matrix);
669	clpMatrix->makeSpecialColumnCopy();
670	}
671	}
672	if (method==ClpSolve::useDual) {
673	// dual
674	if (doCrash)
675	solveOptions.setSpecialOption(0,1,doCrash); // crash
676	else if (doIdiot)
677	solveOptions.setSpecialOption(0,2,doIdiot); // possible idiot
678	} else if (method==ClpSolve::usePrimalorSprint) {
679	// primal
680	// if slp turn everything off
681	if (slpValue>0) {
682	doCrash=false;
683	doSprint=0;
684	doIdiot=-1;
685	solveOptions.setSpecialOption(1,10,slpValue); // slp
686	method=ClpSolve::usePrimal;
687	}
688	if (doCrash) {
689	solveOptions.setSpecialOption(1,1,doCrash); // crash
690	} else if (doSprint>0) {
691	// sprint overrides idiot
692	solveOptions.setSpecialOption(1,3,doSprint); // sprint
693	} else if (doIdiot>0) {
694	solveOptions.setSpecialOption(1,2,doIdiot); // idiot
695	} else if (slpValue<=0) {
696	if (doIdiot==0) {
697	if (doSprint==0)
698	solveOptions.setSpecialOption(1,4); // all slack
699	else
700	solveOptions.setSpecialOption(1,9); // all slack or sprint
701	} else {
702	if (doSprint==0)
703	solveOptions.setSpecialOption(1,8); // all slack or idiot
704	else
705	solveOptions.setSpecialOption(1,7); // initiative
706	}
707	}
708	if (basisHasValues==-1)
709	solveOptions.setSpecialOption(1,11); // switch off values
710	} else if (method==ClpSolve::useBarrier\|\|method==ClpSolve::useBarrierNoCross) {
711	int barrierOptions = choleskyType;
712	if (scaleBarrier)
713	barrierOptions \|= 8;
714	if (doKKT)
715	barrierOptions \|= 16;
716	if (gamma)
717	barrierOptions \|= 32*gamma;
718	if (crossover==3)
719	barrierOptions \|= 256; // try presolve in crossover
720	solveOptions.setSpecialOption(4,barrierOptions);
721	}
722	int status;
723	if (cppValue>=0) {
724	// generate code
725	FILE * fp = fopen("user_driver.cpp","w");
726	if (fp) {
727	// generate enough to do solveOptions
728	model2->generateCpp(fp);
729	solveOptions.generateCpp(fp);
730	fclose(fp);
731	// now call generate code
732	generateCode("user_driver.cpp",cppValue);
733	} else {
734	std::cout<<"Unable to open file user_driver.cpp"<<std::endl;
735	}
736	}
737	#ifdef CLP_MULTIPLE_FACTORIZATIONS
738	int denseCode = parameters[whichParam(DENSE,numberParameters,parameters)].intValue();
739	model2->factorization()->setGoDenseThreshold(denseCode);
740	int smallCode = parameters[whichParam(SMALLFACT,numberParameters,parameters)].intValue();
741	model2->factorization()->setGoSmallThreshold(smallCode);
742	model2->factorization()->goDenseOrSmall(model2->numberRows());
743	#endif
744	try {
745	status=model2->initialSolve(solveOptions);
746	}
747	catch (CoinError e) {
748	e.print();
749	status=-1;
750	}
751	if (dualize) {
752	int returnCode=static_cast<ClpSimplexOther *> (models+iModel)->restoreFromDual(model2);
753	if (model2->status()==3)
754	returnCode=0;
755	delete model2;
756	if (returnCode&&dualize!=2) {
757	currentModel = models+iModel;
758	// register signal handler
759	signal(SIGINT,signal_handler);
760	models[iModel].primal(1);
761	currentModel=NULL;
762	}
763	}
764	if (status>=0)
765	basisHasValues=1;
766	} else {
767	std::cout<<"** Current model not valid"<<std::endl;
768	}
769	break;
770	case STATISTICS:
771	if (goodModels[iModel]) {
772	// If presolve on look at presolved
773	bool deleteModel2=false;
774	ClpSimplex * model2 = models+iModel;
775	if (preSolve) {
776	ClpPresolve pinfo;
777	int presolveOptions2 = presolveOptions&~0x40000000;
778	if ((presolveOptions2&0xffff)!=0)
779	pinfo.setPresolveActions(presolveOptions2);
780	pinfo.setSubstitution(substitution);
781	if ((printOptions&1)!=0)
782	pinfo.statistics();
783	double presolveTolerance =
784	parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
785	model2 =
786	pinfo.presolvedModel(models[iModel],presolveTolerance,
787	true,preSolve);
788	if (model2) {
789	printf("Statistics for presolved model\n");
790	deleteModel2=true;
791	} else {
792	printf("Presolved model looks infeasible - will use unpresolved\n");
793	model2 = models+iModel;
794	}
795	} else {
796	printf("Statistics for unpresolved model\n");
797	model2 = models+iModel;
798	}
799	statistics(models+iModel,model2);
800	if (deleteModel2)
801	delete model2;
802	} else {
803	std::cout<<"** Current model not valid"<<std::endl;
804	}
805	break;
806	case TIGHTEN:
807	if (goodModels[iModel]) {
808	int numberInfeasibilities = models[iModel].tightenPrimalBounds();
809	if (numberInfeasibilities)
810	std::cout<<"** Analysis indicates model infeasible"<<std::endl;
811	} else {
812	std::cout<<"** Current model not valid"<<std::endl;
813	}
814	break;
815	case PLUSMINUS:
816	if (goodModels[iModel]) {
817	ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
818	ClpPackedMatrix* clpMatrix =
819	dynamic_cast< ClpPackedMatrix*>(saveMatrix);
820	if (clpMatrix) {
821	ClpPlusMinusOneMatrix * newMatrix = new ClpPlusMinusOneMatrix(*(clpMatrix->matrix()));
822	if (newMatrix->getIndices()) {
823	models[iModel].replaceMatrix(newMatrix);
824	delete saveMatrix;
825	std::cout<<"Matrix converted to +- one matrix"<<std::endl;
826	} else {
827	std::cout<<"Matrix can not be converted to +- 1 matrix"<<std::endl;
828	}
829	} else {
830	std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
831	}
832	} else {
833	std::cout<<"** Current model not valid"<<std::endl;
834	}
835	break;
836	case NETWORK:
837	if (goodModels[iModel]) {
838	ClpMatrixBase * saveMatrix = models[iModel].clpMatrix();
839	ClpPackedMatrix* clpMatrix =
840	dynamic_cast< ClpPackedMatrix*>(saveMatrix);
841	if (clpMatrix) {
842	ClpNetworkMatrix * newMatrix = new ClpNetworkMatrix(*(clpMatrix->matrix()));
843	if (newMatrix->getIndices()) {
844	models[iModel].replaceMatrix(newMatrix);
845	delete saveMatrix;
846	std::cout<<"Matrix converted to network matrix"<<std::endl;
847	} else {
848	std::cout<<"Matrix can not be converted to network matrix"<<std::endl;
849	}
850	} else {
851	std::cout<<"Matrix not a ClpPackedMatrix"<<std::endl;
852	}
853	} else {
854	std::cout<<"** Current model not valid"<<std::endl;
855	}
856	break;
857	case IMPORT:
858	{
859	// get next field
860	field = CoinReadGetString(argc,argv);
861	if (field=="$") {
862	field = parameters[iParam].stringValue();
863	} else if (field=="EOL") {
864	parameters[iParam].printString();
865	break;
866	} else {
867	parameters[iParam].setStringValue(field);
868	}
869	std::string fileName;
870	bool canOpen=false;
871	// See if gmpl file
872	int gmpl=0;
873	std::string gmplData;
874	if (field=="-") {
875	// stdin
876	canOpen=true;
877	fileName = "-";
878	} else {
879	// See if .lp
880	{
881	const char * c_name = field.c_str();
882	int length = strlen(c_name);
883	if (length>3&&!strncmp(c_name+length-3,".lp",3))
884	gmpl=-1; // .lp
885	}
886	bool absolutePath;
887	if (dirsep=='/') {
888	// non Windows (or cygwin)
889	absolutePath=(field[0]=='/');
890	} else {
891	//Windows (non cycgwin)
892	absolutePath=(field[0]=='\\');
893	// but allow for :
894	if (strchr(field.c_str(),':'))
895	absolutePath=true;
896	}
897	if (absolutePath) {
898	fileName = field;
899	} else if (field[0]=='~') {
900	char * environVar = getenv("HOME");
901	if (environVar) {
902	std::string home(environVar);
903	field=field.erase(0,1);
904	fileName = home+field;
905	} else {
906	fileName=field;
907	}
908	} else {
909	fileName = directory+field;
910	// See if gmpl (model & data) - or even lp file
911	int length = field.size();
912	int percent = field.find('%');
913	if (percent<length&&percent>0) {
914	gmpl=1;
915	fileName = directory+field.substr(0,percent);
916	gmplData = directory+field.substr(percent+1);
917	if (percent<length-1)
918	gmpl=2; // two files
919	printf("GMPL model file %s and data file %s\n",
920	fileName.c_str(),gmplData.c_str());
921	}
922	}
923	std::string name=fileName;
924	if (fileCoinReadable(name)) {
925	// can open - lets go for it
926	canOpen=true;
927	if (gmpl==2) {
928	FILE *fp;
929	fp=fopen(gmplData.c_str(),"r");
930	if (fp) {
931	fclose(fp);
932	} else {
933	canOpen=false;
934	std::cout<<"Unable to open file "<<gmplData<<std::endl;
935	}
936	}
937	} else {
938	std::cout<<"Unable to open file "<<fileName<<std::endl;
939	}
940	}
941	if (canOpen) {
942	int status;
943	if (!gmpl)
944	status =models[iModel].readMps(fileName.c_str(),
945	keepImportNames!=0,
946	allowImportErrors!=0);
947	else if (gmpl>0)
948	status= models[iModel].readGMPL(fileName.c_str(),
949	(gmpl==2) ? gmplData.c_str() : NULL,
950	keepImportNames!=0);
951	else
952	status= models[iModel].readLp(fileName.c_str(),1.0e-12);
953	if (!status\|\|(status>0&&allowImportErrors)) {
954	goodModels[iModel]=true;
955	// sets to all slack (not necessary?)
956	models[iModel].createStatus();
957	time2 = CoinCpuTime();
958	totalTime += time2-time1;
959	time1=time2;
960	// Go to canned file if just input file
961	if (CbcOrClpRead_mode==2&&argc==2) {
962	// only if ends .mps
963	char * find = const_cast<char *>(strstr(fileName.c_str(),".mps"));
964	if (find&&find[4]=='\0') {
965	find[1]='p'; find[2]='a';find[3]='r';
966	FILE *fp=fopen(fileName.c_str(),"r");
967	if (fp) {
968	CbcOrClpReadCommand=fp; // Read from that file
969	CbcOrClpRead_mode=-1;
970	}
971	}
972	}
973	} else {
974	// errors
975	std::cout<<"There were "<<status<<
976	" errors on input"<<std::endl;
977	}
978	}
979	}
980	break;
981	case EXPORT:
982	if (goodModels[iModel]) {
983	double objScale =
984	parameters[whichParam(OBJSCALE2,numberParameters,parameters)].doubleValue();
985	if (objScale!=1.0) {
986	int iColumn;
987	int numberColumns=models[iModel].numberColumns();
988	double * dualColumnSolution =
989	models[iModel].dualColumnSolution();
990	ClpObjective * obj = models[iModel].objectiveAsObject();
991	assert(dynamic_cast<ClpLinearObjective *> (obj));
992	double offset;
993	double * objective = obj->gradient(NULL,NULL,offset,true);
994	for (iColumn=0;iColumn<numberColumns;iColumn++) {
995	dualColumnSolution[iColumn] *= objScale;
996	objective[iColumn] *= objScale;;
997	}
998	int iRow;
999	int numberRows=models[iModel].numberRows();
1000	double * dualRowSolution =
1001	models[iModel].dualRowSolution();
1002	for (iRow=0;iRow<numberRows;iRow++)
1003	dualRowSolution[iRow] *= objScale;
1004	models[iModel].setObjectiveOffset(objScale*models[iModel].objectiveOffset());
1005	}
1006	// get next field
1007	field = CoinReadGetString(argc,argv);
1008	if (field=="$") {
1009	field = parameters[iParam].stringValue();
1010	} else if (field=="EOL") {
1011	parameters[iParam].printString();
1012	break;
1013	} else {
1014	parameters[iParam].setStringValue(field);
1015	}
1016	std::string fileName;
1017	bool canOpen=false;
1018	if (field[0]=='/'\|\|field[0]=='\\') {
1019	fileName = field;
1020	} else if (field[0]=='~') {
1021	char * environVar = getenv("HOME");
1022	if (environVar) {
1023	std::string home(environVar);
1024	field=field.erase(0,1);
1025	fileName = home+field;
1026	} else {
1027	fileName=field;
1028	}
1029	} else {
1030	fileName = directory+field;
1031	}
1032	FILE *fp=fopen(fileName.c_str(),"w");
1033	if (fp) {
1034	// can open - lets go for it
1035	fclose(fp);
1036	canOpen=true;
1037	} else {
1038	std::cout<<"Unable to open file "<<fileName<<std::endl;
1039	}
1040	if (canOpen) {
1041	// If presolve on then save presolved
1042	bool deleteModel2=false;
1043	ClpSimplex * model2 = models+iModel;
1044	if (dualize&&dualize<3) {
1045	model2 = static_cast<ClpSimplexOther *> (model2)->dualOfModel();
1046	printf("Dual of model has %d rows and %d columns\n",
1047	model2->numberRows(),model2->numberColumns());
1048	model2->setOptimizationDirection(1.0);
1049	preSolve=0; // as picks up from model
1050	}
1051	if (preSolve) {
1052	ClpPresolve pinfo;
1053	int presolveOptions2 = presolveOptions&~0x40000000;
1054	if ((presolveOptions2&0xffff)!=0)
1055	pinfo.setPresolveActions(presolveOptions2);
1056	pinfo.setSubstitution(substitution);
1057	if ((printOptions&1)!=0)
1058	pinfo.statistics();
1059	double presolveTolerance =
1060	parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
1061	model2 =
1062	pinfo.presolvedModel(models[iModel],presolveTolerance,
1063	true,preSolve,false,false);
1064	if (model2) {
1065	printf("Saving presolved model on %s\n",
1066	fileName.c_str());
1067	deleteModel2=true;
1068	} else {
1069	printf("Presolved model looks infeasible - saving original on %s\n",
1070	fileName.c_str());
1071	deleteModel2=false;
1072	model2 = models+iModel;
1073
1074	}
1075	} else {
1076	printf("Saving model on %s\n",
1077	fileName.c_str());
1078	}
1079	#if 0
1080	// Convert names
1081	int iRow;
1082	int numberRows=model2->numberRows();
1083	int iColumn;
1084	int numberColumns=model2->numberColumns();
1085
1086	char ** rowNames = NULL;
1087	char ** columnNames = NULL;
1088	if (model2->lengthNames()) {
1089	rowNames = new char * [numberRows];
1090	for (iRow=0;iRow<numberRows;iRow++) {
1091	rowNames[iRow] =
1092	CoinStrdup(model2->rowName(iRow).c_str());
1093	#ifdef STRIPBLANKS
1094	char * xx = rowNames[iRow];
1095	int i;
1096	int length = strlen(xx);
1097	int n=0;
1098	for (i=0;i<length;i++) {
1099	if (xx[i]!=' ')
1100	xx[n++]=xx[i];
1101	}
1102	xx[n]='\0';
1103	#endif
1104	}
1105
1106	columnNames = new char * [numberColumns];
1107	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1108	columnNames[iColumn] =
1109	CoinStrdup(model2->columnName(iColumn).c_str());
1110	#ifdef STRIPBLANKS
1111	char * xx = columnNames[iColumn];
1112	int i;
1113	int length = strlen(xx);
1114	int n=0;
1115	for (i=0;i<length;i++) {
1116	if (xx[i]!=' ')
1117	xx[n++]=xx[i];
1118	}
1119	xx[n]='\0';
1120	#endif
1121	}
1122	}
1123	CoinMpsIO writer;
1124	writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
1125	model2->getColLower(), model2->getColUpper(),
1126	model2->getObjCoefficients(),
1127	(const char) 0 /integrality*/,
1128	model2->getRowLower(), model2->getRowUpper(),
1129	columnNames, rowNames);
1130	// Pass in array saying if each variable integer
1131	writer.copyInIntegerInformation(model2->integerInformation());
1132	writer.setObjectiveOffset(model2->objectiveOffset());
1133	writer.writeMps(fileName.c_str(),0,1,1);
1134	if (rowNames) {
1135	for (iRow=0;iRow<numberRows;iRow++) {
1136	free(rowNames[iRow]);
1137	}
1138	delete [] rowNames;
1139	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1140	free(columnNames[iColumn]);
1141	}
1142	delete [] columnNames;
1143	}
1144	#else
1145	model2->writeMps(fileName.c_str(),(outputFormat-1)/2,1+((outputFormat-1)&1));
1146	#endif
1147	if (deleteModel2)
1148	delete model2;
1149	time2 = CoinCpuTime();
1150	totalTime += time2-time1;
1151	time1=time2;
1152	}
1153	} else {
1154	std::cout<<"** Current model not valid"<<std::endl;
1155	}
1156	break;
1157	case BASISIN:
1158	if (goodModels[iModel]) {
1159	// get next field
1160	field = CoinReadGetString(argc,argv);
1161	if (field=="$") {
1162	field = parameters[iParam].stringValue();
1163	} else if (field=="EOL") {
1164	parameters[iParam].printString();
1165	break;
1166	} else {
1167	parameters[iParam].setStringValue(field);
1168	}
1169	std::string fileName;
1170	bool canOpen=false;
1171	if (field=="-") {
1172	// stdin
1173	canOpen=true;
1174	fileName = "-";
1175	} else {
1176	if (field[0]=='/'\|\|field[0]=='\\') {
1177	fileName = field;
1178	} else if (field[0]=='~') {
1179	char * environVar = getenv("HOME");
1180	if (environVar) {
1181	std::string home(environVar);
1182	field=field.erase(0,1);
1183	fileName = home+field;
1184	} else {
1185	fileName=field;
1186	}
1187	} else {
1188	fileName = directory+field;
1189	}
1190	FILE *fp=fopen(fileName.c_str(),"r");
1191	if (fp) {
1192	// can open - lets go for it
1193	fclose(fp);
1194	canOpen=true;
1195	} else {
1196	std::cout<<"Unable to open file "<<fileName<<std::endl;
1197	}
1198	}
1199	if (canOpen) {
1200	int values = models[iModel].readBasis(fileName.c_str());
1201	if (values==0)
1202	basisHasValues=-1;
1203	else
1204	basisHasValues=1;
1205	}
1206	} else {
1207	std::cout<<"** Current model not valid"<<std::endl;
1208	}
1209	break;
1210	case PRINTMASK:
1211	// get next field
1212	{
1213	std::string name = CoinReadGetString(argc,argv);
1214	if (name!="EOL") {
1215	parameters[iParam].setStringValue(name);
1216	printMask = name;
1217	} else {
1218	parameters[iParam].printString();
1219	}
1220	}
1221	break;
1222	case BASISOUT:
1223	if (goodModels[iModel]) {
1224	// get next field
1225	field = CoinReadGetString(argc,argv);
1226	if (field=="$") {
1227	field = parameters[iParam].stringValue();
1228	} else if (field=="EOL") {
1229	parameters[iParam].printString();
1230	break;
1231	} else {
1232	parameters[iParam].setStringValue(field);
1233	}
1234	std::string fileName;
1235	bool canOpen=false;
1236	if (field[0]=='/'\|\|field[0]=='\\') {
1237	fileName = field;
1238	} else if (field[0]=='~') {
1239	char * environVar = getenv("HOME");
1240	if (environVar) {
1241	std::string home(environVar);
1242	field=field.erase(0,1);
1243	fileName = home+field;
1244	} else {
1245	fileName=field;
1246	}
1247	} else {
1248	fileName = directory+field;
1249	}
1250	FILE *fp=fopen(fileName.c_str(),"w");
1251	if (fp) {
1252	// can open - lets go for it
1253	fclose(fp);
1254	canOpen=true;
1255	} else {
1256	std::cout<<"Unable to open file "<<fileName<<std::endl;
1257	}
1258	if (canOpen) {
1259	ClpSimplex * model2 = models+iModel;
1260	model2->writeBasis(fileName.c_str(),outputFormat>1,outputFormat-2);
1261	time2 = CoinCpuTime();
1262	totalTime += time2-time1;
1263	time1=time2;
1264	}
1265	} else {
1266	std::cout<<"** Current model not valid"<<std::endl;
1267	}
1268	break;
1269	case SAVE:
1270	{
1271	// get next field
1272	field = CoinReadGetString(argc,argv);
1273	if (field=="$") {
1274	field = parameters[iParam].stringValue();
1275	} else if (field=="EOL") {
1276	parameters[iParam].printString();
1277	break;
1278	} else {
1279	parameters[iParam].setStringValue(field);
1280	}
1281	std::string fileName;
1282	bool canOpen=false;
1283	if (field[0]=='/'\|\|field[0]=='\\') {
1284	fileName = field;
1285	} else if (field[0]=='~') {
1286	char * environVar = getenv("HOME");
1287	if (environVar) {
1288	std::string home(environVar);
1289	field=field.erase(0,1);
1290	fileName = home+field;
1291	} else {
1292	fileName=field;
1293	}
1294	} else {
1295	fileName = directory+field;
1296	}
1297	FILE *fp=fopen(fileName.c_str(),"wb");
1298	if (fp) {
1299	// can open - lets go for it
1300	fclose(fp);
1301	canOpen=true;
1302	} else {
1303	std::cout<<"Unable to open file "<<fileName<<std::endl;
1304	}
1305	if (canOpen) {
1306	int status;
1307	// If presolve on then save presolved
1308	bool deleteModel2=false;
1309	ClpSimplex * model2 = models+iModel;
1310	if (preSolve) {
1311	ClpPresolve pinfo;
1312	double presolveTolerance =
1313	parameters[whichParam(PRESOLVETOLERANCE,numberParameters,parameters)].doubleValue();
1314	model2 =
1315	pinfo.presolvedModel(models[iModel],presolveTolerance,
1316	false,preSolve);
1317	if (model2) {
1318	printf("Saving presolved model on %s\n",
1319	fileName.c_str());
1320	deleteModel2=true;
1321	} else {
1322	printf("Presolved model looks infeasible - saving original on %s\n",
1323	fileName.c_str());
1324	deleteModel2=false;
1325	model2 = models+iModel;
1326
1327	}
1328	} else {
1329	printf("Saving model on %s\n",
1330	fileName.c_str());
1331	}
1332	status =model2->saveModel(fileName.c_str());
1333	if (deleteModel2)
1334	delete model2;
1335	if (!status) {
1336	goodModels[iModel]=true;
1337	time2 = CoinCpuTime();
1338	totalTime += time2-time1;
1339	time1=time2;
1340	} else {
1341	// errors
1342	std::cout<<"There were errors on output"<<std::endl;
1343	}
1344	}
1345	}
1346	break;
1347	case RESTORE:
1348	{
1349	// get next field
1350	field = CoinReadGetString(argc,argv);
1351	if (field=="$") {
1352	field = parameters[iParam].stringValue();
1353	} else if (field=="EOL") {
1354	parameters[iParam].printString();
1355	break;
1356	} else {
1357	parameters[iParam].setStringValue(field);
1358	}
1359	std::string fileName;
1360	bool canOpen=false;
1361	if (field[0]=='/'\|\|field[0]=='\\') {
1362	fileName = field;
1363	} else if (field[0]=='~') {
1364	char * environVar = getenv("HOME");
1365	if (environVar) {
1366	std::string home(environVar);
1367	field=field.erase(0,1);
1368	fileName = home+field;
1369	} else {
1370	fileName=field;
1371	}
1372	} else {
1373	fileName = directory+field;
1374	}
1375	FILE *fp=fopen(fileName.c_str(),"rb");
1376	if (fp) {
1377	// can open - lets go for it
1378	fclose(fp);
1379	canOpen=true;
1380	} else {
1381	std::cout<<"Unable to open file "<<fileName<<std::endl;
1382	}
1383	if (canOpen) {
1384	int status =models[iModel].restoreModel(fileName.c_str());
1385	if (!status) {
1386	goodModels[iModel]=true;
1387	time2 = CoinCpuTime();
1388	totalTime += time2-time1;
1389	time1=time2;
1390	} else {
1391	// errors
1392	std::cout<<"There were errors on input"<<std::endl;
1393	}
1394	}
1395	}
1396	break;
1397	case MAXIMIZE:
1398	models[iModel].setOptimizationDirection(-1);
1399	break;
1400	case MINIMIZE:
1401	models[iModel].setOptimizationDirection(1);
1402	break;
1403	case ALLSLACK:
1404	models[iModel].allSlackBasis(true);
1405	break;
1406	case REVERSE:
1407	if (goodModels[iModel]) {
1408	int iColumn;
1409	int numberColumns=models[iModel].numberColumns();
1410	double * dualColumnSolution =
1411	models[iModel].dualColumnSolution();
1412	ClpObjective * obj = models[iModel].objectiveAsObject();
1413	assert(dynamic_cast<ClpLinearObjective *> (obj));
1414	double offset;
1415	double * objective = obj->gradient(NULL,NULL,offset,true);
1416	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1417	dualColumnSolution[iColumn] = -dualColumnSolution[iColumn];
1418	objective[iColumn] = -objective[iColumn];
1419	}
1420	int iRow;
1421	int numberRows=models[iModel].numberRows();
1422	double * dualRowSolution =
1423	models[iModel].dualRowSolution();
1424	for (iRow=0;iRow<numberRows;iRow++) {
1425	dualRowSolution[iRow] = -dualRowSolution[iRow];
1426	}
1427	models[iModel].setObjectiveOffset(-models[iModel].objectiveOffset());
1428	}
1429	break;
1430	case DIRECTORY:
1431	{
1432	std::string name = CoinReadGetString(argc,argv);
1433	if (name!="EOL") {
1434	int length=name.length();
1435	if (name[length-1]==dirsep) {
1436	directory = name;
1437	} else {
1438	directory = name+dirsep;
1439	}
1440	parameters[iParam].setStringValue(directory);
1441	} else {
1442	parameters[iParam].printString();
1443	}
1444	}
1445	break;
1446	case DIRSAMPLE:
1447	{
1448	std::string name = CoinReadGetString(argc,argv);
1449	if (name!="EOL") {
1450	int length=name.length();
1451	if (name[length-1]==dirsep) {
1452	dirSample = name;
1453	} else {
1454	dirSample = name+dirsep;
1455	}
1456	parameters[iParam].setStringValue(dirSample);
1457	} else {
1458	parameters[iParam].printString();
1459	}
1460	}
1461	break;
1462	case DIRNETLIB:
1463	{
1464	std::string name = CoinReadGetString(argc,argv);
1465	if (name!="EOL") {
1466	int length=name.length();
1467	if (name[length-1]==dirsep) {
1468	dirNetlib = name;
1469	} else {
1470	dirNetlib = name+dirsep;
1471	}
1472	parameters[iParam].setStringValue(dirNetlib);
1473	} else {
1474	parameters[iParam].printString();
1475	}
1476	}
1477	break;
1478	case DIRMIPLIB:
1479	{
1480	std::string name = CoinReadGetString(argc,argv);
1481	if (name!="EOL") {
1482	int length=name.length();
1483	if (name[length-1]==dirsep) {
1484	dirMiplib = name;
1485	} else {
1486	dirMiplib = name+dirsep;
1487	}
1488	parameters[iParam].setStringValue(dirMiplib);
1489	} else {
1490	parameters[iParam].printString();
1491	}
1492	}
1493	break;
1494	case STDIN:
1495	CbcOrClpRead_mode=-1;
1496	break;
1497	case NETLIB_DUAL:
1498	case NETLIB_EITHER:
1499	case NETLIB_BARRIER:
1500	case NETLIB_PRIMAL:
1501	case NETLIB_TUNE:
1502	{
1503	// create fields for unitTest
1504	const char * fields[4];
1505	int nFields=4;
1506	fields[0]="fake main from unitTest";
1507	std::string mpsfield = "-dirSample=";
1508	mpsfield += dirSample.c_str();
1509	fields[1]=mpsfield.c_str();
1510	std::string netfield = "-dirNetlib=";
1511	netfield += dirNetlib.c_str();
1512	fields[2]=netfield.c_str();
1513	fields[3]="-netlib";
1514	int algorithm;
1515	if (type==NETLIB_DUAL) {
1516	std::cerr<<"Doing netlib with dual algorithm"<<std::endl;
1517	algorithm =0;
1518	} else if (type==NETLIB_BARRIER) {
1519	std::cerr<<"Doing netlib with barrier algorithm"<<std::endl;
1520	algorithm =2;
1521	} else if (type==NETLIB_EITHER) {
1522	std::cerr<<"Doing netlib with dual or primal algorithm"<<std::endl;
1523	algorithm =3;
1524	} else if (type==NETLIB_TUNE) {
1525	std::cerr<<"Doing netlib with best algorithm!"<<std::endl;
1526	algorithm =5;
1527	// uncomment next to get active tuning
1528	// algorithm=6;
1529	} else {
1530	std::cerr<<"Doing netlib with primal agorithm"<<std::endl;
1531	algorithm=1;
1532	}
1533	int specialOptions = models[iModel].specialOptions();
1534	models[iModel].setSpecialOptions(0);
1535	mainTest(nFields,fields,algorithm,models[iModel],
1536	(preSolve!=0),specialOptions,doVector!=0);
1537	}
1538	break;
1539	case UNITTEST:
1540	{
1541	// create fields for unitTest
1542	const char * fields[2];
1543	int nFields=2;
1544	fields[0]="fake main from unitTest";
1545	std::string dirfield = "-dirSample=";
1546	dirfield += dirSample.c_str();
1547	fields[1]=dirfield.c_str();
1548	int specialOptions = models[iModel].specialOptions();
1549	models[iModel].setSpecialOptions(0);
1550	int algorithm=-1;
1551	if (models[iModel].numberRows())
1552	algorithm=7;
1553	mainTest(nFields,fields,algorithm,models[iModel],(preSolve!=0),specialOptions,doVector!=0);
1554	}
1555	break;
1556	case FAKEBOUND:
1557	if (goodModels[iModel]) {
1558	// get bound
1559	double value = CoinReadGetDoubleField(argc,argv,&valid);
1560	if (!valid) {
1561	std::cout<<"Setting "<<parameters[iParam].name()<<
1562	" to DEBUG "<<value<<std::endl;
1563	int iRow;
1564	int numberRows=models[iModel].numberRows();
1565	double * rowLower = models[iModel].rowLower();
1566	double * rowUpper = models[iModel].rowUpper();
1567	for (iRow=0;iRow<numberRows;iRow++) {
1568	// leave free ones for now
1569	if (rowLower[iRow]>-1.0e20\|\|rowUpper[iRow]<1.0e20) {
1570	rowLower[iRow]=CoinMax(rowLower[iRow],-value);
1571	rowUpper[iRow]=CoinMin(rowUpper[iRow],value);
1572	}
1573	}
1574	int iColumn;
1575	int numberColumns=models[iModel].numberColumns();
1576	double * columnLower = models[iModel].columnLower();
1577	double * columnUpper = models[iModel].columnUpper();
1578	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1579	// leave free ones for now
1580	if (columnLower[iColumn]>-1.0e20\|\|
1581	columnUpper[iColumn]<1.0e20) {
1582	columnLower[iColumn]=CoinMax(columnLower[iColumn],-value);
1583	columnUpper[iColumn]=CoinMin(columnUpper[iColumn],value);
1584	}
1585	}
1586	} else if (valid==1) {
1587	abort();
1588	} else {
1589	std::cout<<"enter value for "<<parameters[iParam].name()<<
1590	std::endl;
1591	}
1592	}
1593	break;
1594	case REALLY_SCALE:
1595	if (goodModels[iModel]) {
1596	ClpSimplex newModel(models[iModel],
1597	models[iModel].scalingFlag());
1598	printf("model really really scaled\n");
1599	models[iModel]=newModel;
1600	}
1601	break;
1602	case USERCLP:
1603	// Replace the sample code by whatever you want
1604	if (goodModels[iModel]) {
1605	ClpSimplex * thisModel = &models[iModel];
1606	printf("Dummy user code - model has %d rows and %d columns\n",
1607	thisModel->numberRows(),thisModel->numberColumns());
1608	}
1609	break;
1610	case HELP:
1611	std::cout<<"Coin LP version "<<CLPVERSION
1612	<<", build "<<__DATE__<<std::endl;
1613	std::cout<<"Non default values:-"<<std::endl;
1614	std::cout<<"Perturbation "<<models[0].perturbation()<<" (default 100)"
1615	<<std::endl;
1616	CoinReadPrintit(
1617	"Presolve being done with 5 passes\n\
1618	Dual steepest edge steep/partial on matrix shape and factorization density\n\
1619	Clpnnnn taken out of messages\n\
1620	If Factorization frequency default then done on size of matrix\n\n\
1621	(-)unitTest, (-)netlib or (-)netlibp will do standard tests\n\n\
1622	You can switch to interactive mode at any time so\n\
1623	clp watson.mps -scaling off -primalsimplex\nis the same as\n\
1624	clp watson.mps -\nscaling off\nprimalsimplex"
1625	);
1626	break;
1627	case SOLUTION:
1628	if (goodModels[iModel]) {
1629	// get next field
1630	field = CoinReadGetString(argc,argv);
1631	if (field=="$") {
1632	field = parameters[iParam].stringValue();
1633	} else if (field=="EOL") {
1634	parameters[iParam].printString();
1635	break;
1636	} else {
1637	parameters[iParam].setStringValue(field);
1638	}
1639	std::string fileName;
1640	FILE *fp=NULL;
1641	if (field=="-"\|\|field=="EOL"\|\|field=="stdout") {
1642	// stdout
1643	fp=stdout;
1644	fprintf(fp,"\n");
1645	} else if (field=="stderr") {
1646	// stderr
1647	fp=stderr;
1648	fprintf(fp,"\n");
1649	} else {
1650	if (field[0]=='/'\|\|field[0]=='\\') {
1651	fileName = field;
1652	} else if (field[0]=='~') {
1653	char * environVar = getenv("HOME");
1654	if (environVar) {
1655	std::string home(environVar);
1656	field=field.erase(0,1);
1657	fileName = home+field;
1658	} else {
1659	fileName=field;
1660	}
1661	} else {
1662	fileName = directory+field;
1663	}
1664	fp=fopen(fileName.c_str(),"w");
1665	}
1666	if (fp) {
1667	// Write solution header (suggested by Luigi Poderico)
1668	double objValue = models[iModel].getObjValue()*models[iModel].getObjSense();
1669	int iStat = models[iModel].status();
1670	if (iStat==0) {
1671	fprintf(fp, "optimal\n" );
1672	} else if (iStat==1) {
1673	// infeasible
1674	fprintf(fp, "infeasible\n" );
1675	} else if (iStat==2) {
1676	// unbounded
1677	fprintf(fp, "unbounded\n" );
1678	} else if (iStat==3) {
1679	fprintf(fp, "stopped on iterations or time\n" );
1680	} else if (iStat==4) {
1681	fprintf(fp, "stopped on difficulties\n" );
1682	} else if (iStat==5) {
1683	fprintf(fp, "stopped on ctrl-c\n" );
1684	} else {
1685	fprintf(fp, "status unknown\n" );
1686	}
1687	fprintf(fp, "Objective value %15.8g\n", objValue);
1688	// make fancy later on
1689	int iRow;
1690	int numberRows=models[iModel].numberRows();
1691	int lengthName = models[iModel].lengthNames(); // 0 if no names
1692	// in general I don't want to pass around massive
1693	// amounts of data but seems simpler here
1694	std::vector<std::string> rowNames =
1695	*(models[iModel].rowNames());
1696	std::vector<std::string> columnNames =
1697	*(models[iModel].columnNames());
1698
1699	double * dualRowSolution = models[iModel].dualRowSolution();
1700	double * primalRowSolution =
1701	models[iModel].primalRowSolution();
1702	double * rowLower = models[iModel].rowLower();
1703	double * rowUpper = models[iModel].rowUpper();
1704	double primalTolerance = models[iModel].primalTolerance();
1705	char format[6];
1706	sprintf(format,"%%-%ds",CoinMax(lengthName,8));
1707	bool doMask = (printMask!=""&&lengthName);
1708	int * maskStarts=NULL;
1709	int maxMasks=0;
1710	char ** masks =NULL;
1711	if (doMask) {
1712	int nAst =0;
1713	const char * pMask2 = printMask.c_str();
1714	char pMask[100];
1715	int iChar;
1716	int lengthMask = strlen(pMask2);
1717	assert (lengthMask<100);
1718	if (*pMask2=='"') {
1719	if (pMask2[lengthMask-1]!='"') {
1720	printf("mismatched \" in mask %s\n",pMask2);
1721	break;
1722	} else {
1723	strcpy(pMask,pMask2+1);
1724	*strchr(pMask,'"')='\0';
1725	}
1726	} else if (*pMask2=='\'') {
1727	if (pMask2[lengthMask-1]!='\'') {
1728	printf("mismatched ' in mask %s\n",pMask2);
1729	break;
1730	} else {
1731	strcpy(pMask,pMask2+1);
1732	*strchr(pMask,'\'')='\0';
1733	}
1734	} else {
1735	strcpy(pMask,pMask2);
1736	}
1737	if (lengthMask>lengthName) {
1738	printf("mask %s too long - skipping\n",pMask);
1739	break;
1740	}
1741	maxMasks = 1;
1742	for (iChar=0;iChar<lengthMask;iChar++) {
1743	if (pMask[iChar]=='*') {
1744	nAst++;
1745	maxMasks *= (lengthName+1);
1746	}
1747	}
1748	int nEntries = 1;
1749	maskStarts = new int[lengthName+2];
1750	masks = new char * [maxMasks];
1751	char ** newMasks = new char * [maxMasks];
1752	int i;
1753	for (i=0;i<maxMasks;i++) {
1754	masks[i] = new char[lengthName+1];
1755	newMasks[i] = new char[lengthName+1];
1756	}
1757	strcpy(masks[0],pMask);
1758	for (int iAst=0;iAst<nAst;iAst++) {
1759	int nOldEntries = nEntries;
1760	nEntries=0;
1761	for (int iEntry = 0;iEntry<nOldEntries;iEntry++) {
1762	char * oldMask = masks[iEntry];
1763	char * ast = strchr(oldMask,'*');
1764	assert (ast);
1765	int length = strlen(oldMask)-1;
1766	int nBefore = ast-oldMask;
1767	int nAfter = length-nBefore;
1768	// and add null
1769	nAfter++;
1770	for (int i=0;i<=lengthName-length;i++) {
1771	char * maskOut = newMasks[nEntries];
1772	CoinMemcpyN(oldMask,nBefore,maskOut);
1773	for (int k=0;k<i;k++)
1774	maskOut[k+nBefore]='?';
1775	CoinMemcpyN(ast+1,nAfter,maskOut+nBefore+i);
1776	nEntries++;
1777	assert (nEntries<=maxMasks);
1778	}
1779	}
1780	char ** temp = masks;
1781	masks = newMasks;
1782	newMasks = temp;
1783	}
1784	// Now extend and sort
1785	int * sort = new int[nEntries];
1786	for (i=0;i<nEntries;i++) {
1787	char * maskThis = masks[i];
1788	int length = strlen(maskThis);
1789	while (maskThis[length-1]==' ')
1790	length--;
1791	maskThis[length]='\0';
1792	sort[i]=length;
1793	}
1794	CoinSort_2(sort,sort+nEntries,masks);
1795	int lastLength=-1;
1796	for (i=0;i<nEntries;i++) {
1797	int length = sort[i];
1798	while (length>lastLength)
1799	maskStarts[++lastLength] = i;
1800	}
1801	maskStarts[++lastLength]=nEntries;
1802	delete [] sort;
1803	for (i=0;i<maxMasks;i++)
1804	delete [] newMasks[i];
1805	delete [] newMasks;
1806	}
1807	if (printMode>2) {
1808	for (iRow=0;iRow<numberRows;iRow++) {
1809	int type=printMode-3;
1810	if (primalRowSolution[iRow]>rowUpper[iRow]+primalTolerance\|\|
1811	primalRowSolution[iRow]<rowLower[iRow]-primalTolerance) {
1812	fprintf(fp,"** ");
1813	type=2;
1814	} else if (fabs(primalRowSolution[iRow])>1.0e-8) {
1815	type=1;
1816	} else if (numberRows<50) {
1817	type=3;
1818	}
1819	if (doMask&&!maskMatches(maskStarts,masks,rowNames[iRow]))
1820	type=0;
1821	if (type) {
1822	fprintf(fp,"%7d ",iRow);
1823	if (lengthName)
1824	fprintf(fp,format,rowNames[iRow].c_str());
1825	fprintf(fp,"%15.8g %15.8g\n",primalRowSolution[iRow],
1826	dualRowSolution[iRow]);
1827	}
1828	}
1829	}
1830	int iColumn;
1831	int numberColumns=models[iModel].numberColumns();
1832	double * dualColumnSolution =
1833	models[iModel].dualColumnSolution();
1834	double * primalColumnSolution =
1835	models[iModel].primalColumnSolution();
1836	double * columnLower = models[iModel].columnLower();
1837	double * columnUpper = models[iModel].columnUpper();
1838	for (iColumn=0;iColumn<numberColumns;iColumn++) {
1839	int type=(printMode>3) ? 1 : 0;
1840	if (primalColumnSolution[iColumn]>columnUpper[iColumn]+primalTolerance\|\|
1841	primalColumnSolution[iColumn]<columnLower[iColumn]-primalTolerance) {
1842	fprintf(fp,"** ");
1843	type=2;
1844	} else if (fabs(primalColumnSolution[iColumn])>1.0e-8) {
1845	type=1;
1846	} else if (numberColumns<50) {
1847	type=3;
1848	}
1849	if (doMask&&!maskMatches(maskStarts,masks,
1850	columnNames[iColumn]))
1851	type =0;
1852	if (type) {
1853	fprintf(fp,"%7d ",iColumn);
1854	if (lengthName)
1855	fprintf(fp,format,columnNames[iColumn].c_str());
1856	fprintf(fp,"%15.8g %15.8g\n",
1857	primalColumnSolution[iColumn],
1858	dualColumnSolution[iColumn]);
1859	}
1860	}
1861	if (fp!=stdout)
1862	fclose(fp);
1863	if (masks) {
1864	delete [] maskStarts;
1865	for (int i=0;i<maxMasks;i++)
1866	delete [] masks[i];
1867	delete [] masks;
1868	}
1869	} else {
1870	std::cout<<"Unable to open file "<<fileName<<std::endl;
1871	}
1872	} else {
1873	std::cout<<"** Current model not valid"<<std::endl;
1874
1875	}
1876
1877	break;
1878	case SAVESOL:
1879	if (goodModels[iModel]) {
1880	// get next field
1881	field = CoinReadGetString(argc,argv);
1882	if (field=="$") {
1883	field = parameters[iParam].stringValue();
1884	} else if (field=="EOL") {
1885	parameters[iParam].printString();
1886	break;
1887	} else {
1888	parameters[iParam].setStringValue(field);
1889	}
1890	std::string fileName;
1891	if (field[0]=='/'\|\|field[0]=='\\') {
1892	fileName = field;
1893	} else if (field[0]=='~') {
1894	char * environVar = getenv("HOME");
1895	if (environVar) {
1896	std::string home(environVar);
1897	field=field.erase(0,1);
1898	fileName = home+field;
1899	} else {
1900	fileName=field;
1901	}
1902	} else {
1903	fileName = directory+field;
1904	}
1905	saveSolution(models+iModel,fileName);
1906	} else {
1907	std::cout<<"** Current model not valid"<<std::endl;
1908
1909	}
1910	break;
1911	default:
1912	abort();
1913	}
1914	}
1915	} else if (!numberMatches) {
1916	std::cout<<"No match for "<<field<<" - ? for list of commands"
1917	<<std::endl;
1918	} else if (numberMatches==1) {
1919	if (!numberQuery) {
1920	std::cout<<"Short match for "<<field<<" - completion: ";
1921	std::cout<<parameters[firstMatch].matchName()<<std::endl;
1922	} else if (numberQuery) {
1923	std::cout<<parameters[firstMatch].matchName()<<" : ";
1924	std::cout<<parameters[firstMatch].shortHelp()<<std::endl;
1925	if (numberQuery>=2)
1926	parameters[firstMatch].printLongHelp();
1927	}
1928	} else {
1929	if (!numberQuery)
1930	std::cout<<"Multiple matches for "<<field<<" - possible completions:"
1931	<<std::endl;
1932	else
1933	std::cout<<"Completions of "<<field<<":"<<std::endl;
1934	for ( iParam=0; iParam<numberParameters; iParam++ ) {
1935	int match = parameters[iParam].matches(field);
1936	if (match&&parameters[iParam].displayThis()) {
1937	std::cout<<parameters[iParam].matchName();
1938	if (numberQuery>=2)
1939	std::cout<<" : "<<parameters[iParam].shortHelp();
1940	std::cout<<std::endl;
1941	}
1942	}
1943	}
1944	}
1945	delete [] models;
1946	delete [] goodModels;
1947	}
1948	// By now all memory should be freed
1949	#ifdef DMALLOC
1950	dmalloc_log_unfreed();
1951	dmalloc_shutdown();
1952	#endif
1953	return 0;
1954	}
1955	static void breakdown(const char * name, int numberLook, const double * region)
1956	{
1957	double range[] = {
1958	-COIN_DBL_MAX,
1959	-1.0e15,-1.0e11,-1.0e8,-1.0e5,-1.0e4,-1.0e3,-1.0e2,-1.0e1,
1960	-1.0,
1961	-1.0e-1,-1.0e-2,-1.0e-3,-1.0e-4,-1.0e-5,-1.0e-8,-1.0e-11,-1.0e-15,
1962	0.0,
1963	1.0e-15,1.0e-11,1.0e-8,1.0e-5,1.0e-4,1.0e-3,1.0e-2,1.0e-1,
1964	1.0,
1965	1.0e1,1.0e2,1.0e3,1.0e4,1.0e5,1.0e8,1.0e11,1.0e15,
1966	COIN_DBL_MAX};
1967	int nRanges = static_cast<int> (sizeof(range)/sizeof(double));
1968	int * number = new int[nRanges];
1969	memset(number,0,nRanges*sizeof(int));
1970	int * numberExact = new int[nRanges];
1971	memset(numberExact,0,nRanges*sizeof(int));
1972	int i;
1973	for ( i=0;i<numberLook;i++) {
1974	double value = region[i];
1975	for (int j=0;j<nRanges;j++) {
1976	if (value==range[j]) {
1977	numberExact[j]++;
1978	break;
1979	} else if (value<range[j]) {
1980	number[j]++;
1981	break;
1982	}
1983	}
1984	}
1985	printf("\n%s has %d entries\n",name,numberLook);
1986	for (i=0;i<nRanges;i++) {
1987	if (number[i])
1988	printf("%d between %g and %g",number[i],range[i-1],range[i]);
1989	if (numberExact[i]) {
1990	if (number[i])
1991	printf(", ");
1992	printf("%d exactly at %g",numberExact[i],range[i]);
1993	}
1994	if (number[i]+numberExact[i])
1995	printf("\n");
1996	}
1997	delete [] number;
1998	delete [] numberExact;
1999	}
2000	void sortOnOther(int * column,
2001	const CoinBigIndex * rowStart,
2002	int * order,
2003	int * other,
2004	int nRow,
2005	int nInRow,
2006	int where)
2007	{
2008	if (nRow<2\|\|where>=nInRow)
2009	return;
2010	// do initial sort
2011	int kRow;
2012	int iRow;
2013	for ( kRow=0;kRow<nRow;kRow++) {
2014	iRow = order[kRow];
2015	other[kRow]=column[rowStart[iRow]+where];
2016	}
2017	CoinSort_2(other,other+nRow,order);
2018	int first=0;
2019	iRow=order[0];
2020	int firstC=column[rowStart[iRow]+where];
2021	kRow=1;
2022	while (kRow<nRow) {
2023	int lastC=9999999;;
2024	for (;kRow<nRow+1;kRow++) {
2025	if (kRow<nRow) {
2026	iRow=order[kRow];
2027	lastC=column[rowStart[iRow]+where];
2028	} else {
2029	lastC=9999999;
2030	}
2031	if (lastC>firstC)
2032	break;
2033	}
2034	// sort
2035	sortOnOther(column,rowStart,order+first,other,kRow-first,
2036	nInRow,where+1);
2037	firstC=lastC;
2038	first=kRow;
2039	}
2040	}
2041	static void statistics(ClpSimplex * originalModel, ClpSimplex * model)
2042	{
2043	int numberColumns = originalModel->numberColumns();
2044	const char * integerInformation = originalModel->integerInformation();
2045	const double * columnLower = originalModel->columnLower();
2046	const double * columnUpper = originalModel->columnUpper();
2047	int numberIntegers=0;
2048	int numberBinary=0;
2049	int iRow,iColumn;
2050	if (integerInformation) {
2051	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2052	if (integerInformation[iColumn]) {
2053	if (columnUpper[iColumn]>columnLower[iColumn]) {
2054	numberIntegers++;
2055	if (columnUpper[iColumn]==0.0&&columnLower[iColumn]==1)
2056	numberBinary++;
2057	}
2058	}
2059	}
2060	}
2061	numberColumns = model->numberColumns();
2062	int numberRows = model->numberRows();
2063	columnLower = model->columnLower();
2064	columnUpper = model->columnUpper();
2065	const double * rowLower = model->rowLower();
2066	const double * rowUpper = model->rowUpper();
2067	const double * objective = model->objective();
2068	CoinPackedMatrix * matrix = model->matrix();
2069	CoinBigIndex numberElements = matrix->getNumElements();
2070	const int * columnLength = matrix->getVectorLengths();
2071	//const CoinBigIndex * columnStart = matrix->getVectorStarts();
2072	const double * elementByColumn = matrix->getElements();
2073	int * number = new int[numberRows+1];
2074	memset(number,0,(numberRows+1)*sizeof(int));
2075	int numberObjSingletons=0;
2076	/* cType
2077	0 0/inf, 1 0/up, 2 lo/inf, 3 lo/up, 4 free, 5 fix, 6 -inf/0, 7 -inf/up,
2078	8 0/1
2079	*/
2080	int cType[9];
2081	std::string cName[]={"0.0->inf,","0.0->up,","lo->inf,","lo->up,","free,","fixed,","-inf->0.0,",
2082	"-inf->up,","0.0->1.0"};
2083	int nObjective=0;
2084	memset(cType,0,sizeof(cType));
2085	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2086	int length=columnLength[iColumn];
2087	if (length==1&&objective[iColumn])
2088	numberObjSingletons++;
2089	number[length]++;
2090	if (objective[iColumn])
2091	nObjective++;
2092	if (columnLower[iColumn]>-1.0e20) {
2093	if (columnLower[iColumn]==0.0) {
2094	if (columnUpper[iColumn]>1.0e20)
2095	cType[0]++;
2096	else if (columnUpper[iColumn]==1.0)
2097	cType[8]++;
2098	else if (columnUpper[iColumn]==0.0)
2099	cType[5]++;
2100	else
2101	cType[1]++;
2102	} else {
2103	if (columnUpper[iColumn]>1.0e20)
2104	cType[2]++;
2105	else if (columnUpper[iColumn]==columnLower[iColumn])
2106	cType[5]++;
2107	else
2108	cType[3]++;
2109	}
2110	} else {
2111	if (columnUpper[iColumn]>1.0e20)
2112	cType[4]++;
2113	else if (columnUpper[iColumn]==0.0)
2114	cType[6]++;
2115	else
2116	cType[7]++;
2117	}
2118	}
2119	/* rType
2120	0 E 0, 1 E 1, 2 E -1, 3 E other, 4 G 0, 5 G 1, 6 G other,
2121	7 L 0, 8 L 1, 9 L other, 10 Range 0/1, 11 Range other, 12 free
2122	*/
2123	int rType[13];
2124	std::string rName[]={"E 0.0,","E 1.0,","E -1.0,","E other,","G 0.0,","G 1.0,","G other,",
2125	"L 0.0,","L 1.0,","L other,","Range 0.0->1.0,","Range other,","Free"};
2126	memset(rType,0,sizeof(rType));
2127	for (iRow=0;iRow<numberRows;iRow++) {
2128	if (rowLower[iRow]>-1.0e20) {
2129	if (rowLower[iRow]==0.0) {
2130	if (rowUpper[iRow]>1.0e20)
2131	rType[4]++;
2132	else if (rowUpper[iRow]==1.0)
2133	rType[10]++;
2134	else if (rowUpper[iRow]==0.0)
2135	rType[0]++;
2136	else
2137	rType[11]++;
2138	} else if (rowLower[iRow]==1.0) {
2139	if (rowUpper[iRow]>1.0e20)
2140	rType[5]++;
2141	else if (rowUpper[iRow]==rowLower[iRow])
2142	rType[1]++;
2143	else
2144	rType[11]++;
2145	} else if (rowLower[iRow]==-1.0) {
2146	if (rowUpper[iRow]>1.0e20)
2147	rType[6]++;
2148	else if (rowUpper[iRow]==rowLower[iRow])
2149	rType[2]++;
2150	else
2151	rType[11]++;
2152	} else {
2153	if (rowUpper[iRow]>1.0e20)
2154	rType[6]++;
2155	else if (rowUpper[iRow]==rowLower[iRow])
2156	rType[3]++;
2157	else
2158	rType[11]++;
2159	}
2160	} else {
2161	if (rowUpper[iRow]>1.0e20)
2162	rType[12]++;
2163	else if (rowUpper[iRow]==0.0)
2164	rType[7]++;
2165	else if (rowUpper[iRow]==1.0)
2166	rType[8]++;
2167	else
2168	rType[9]++;
2169	}
2170	}
2171	// Basic statistics
2172	printf("\n\nProblem has %d rows, %d columns (%d with objective) and %d elements\n",
2173	numberRows,numberColumns,nObjective,numberElements);
2174	if (number[0]+number[1]) {
2175	printf("There are ");
2176	if (numberObjSingletons)
2177	printf("%d singletons with objective ",numberObjSingletons);
2178	int numberNoObj = number[1]-numberObjSingletons;
2179	if (numberNoObj)
2180	printf("%d singletons with no objective ",numberNoObj);
2181	if (number[0])
2182	printf("** %d columns have no entries",number[0]);
2183	printf("\n");
2184	}
2185	printf("Column breakdown:\n");
2186	int k;
2187	for (k=0;k<static_cast<int> (sizeof(cType)/sizeof(int));k++) {
2188	printf("%d of type %s ",cType[k],cName[k].c_str());
2189	if (((k+1)%3)==0)
2190	printf("\n");
2191	}
2192	if ((k%3)!=0)
2193	printf("\n");
2194	printf("Row breakdown:\n");
2195	for (k=0;k<static_cast<int> (sizeof(rType)/sizeof(int));k++) {
2196	printf("%d of type %s ",rType[k],rName[k].c_str());
2197	if (((k+1)%3)==0)
2198	printf("\n");
2199	}
2200	if ((k%3)!=0)
2201	printf("\n");
2202	#define SYM
2203	#ifndef SYM
2204	if (model->logLevel()<2)
2205	return ;
2206	#endif
2207	int kMax = model->logLevel()>3 ? 1000000 : 10;
2208	k=0;
2209	for (iRow=1;iRow<=numberRows;iRow++) {
2210	if (number[iRow]) {
2211	k++;
2212	printf("%d columns have %d entries\n",number[iRow],iRow);
2213	if (k==kMax)
2214	break;
2215	}
2216	}
2217	if (k<numberRows) {
2218	int kk=k;
2219	k=0;
2220	for (iRow=numberRows;iRow>=1;iRow--) {
2221	if (number[iRow]) {
2222	k++;
2223	if (k==kMax)
2224	break;
2225	}
2226	}
2227	if (k>kk) {
2228	printf("\n .........\n\n");
2229	iRow=k;
2230	k=0;
2231	for (;iRow<numberRows;iRow++) {
2232	if (number[iRow]) {
2233	k++;
2234	printf("%d columns have %d entries\n",number[iRow],iRow);
2235	if (k==kMax)
2236	break;
2237	}
2238	}
2239	}
2240	}
2241	delete [] number;
2242	printf("\n\n");
2243	if (model->logLevel()==63
2244	#ifdef SYM
2245	\|\|true
2246	#endif
2247	) {
2248	// get column copy
2249	CoinPackedMatrix columnCopy = *matrix;
2250	const int * columnLength = columnCopy.getVectorLengths();
2251	number = new int[numberRows+1];
2252	memset(number,0,(numberRows+1)*sizeof(int));
2253	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2254	int length=columnLength[iColumn];
2255	number[length]++;
2256	}
2257	k=0;
2258	for (iRow=1;iRow<=numberRows;iRow++) {
2259	if (number[iRow]) {
2260	k++;
2261	}
2262	}
2263	int * row = columnCopy.getMutableIndices();
2264	const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2265	double * element = columnCopy.getMutableElements();
2266	int * order = new int[numberColumns];
2267	int * other = new int[numberColumns];
2268	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2269	int length=columnLength[iColumn];
2270	order[iColumn]=iColumn;
2271	other[iColumn]=length;
2272	CoinBigIndex start = columnStart[iColumn];
2273	CoinSort_2(row+start,row+start+length,element+start);
2274	}
2275	CoinSort_2(other,other+numberColumns,order);
2276	int jColumn=number[0]+number[1];
2277	for (iRow=2;iRow<=numberRows;iRow++) {
2278	if (number[iRow]) {
2279	printf("XX %d columns have %d entries\n",number[iRow],iRow);
2280	int kColumn = jColumn+number[iRow];
2281	sortOnOther(row,columnStart,
2282	order+jColumn,other,number[iRow],iRow,0);
2283	// Now print etc
2284	if (iRow<500000) {
2285	for (int lColumn =jColumn;lColumn<kColumn;lColumn++) {
2286	iColumn = order[lColumn];
2287	CoinBigIndex start = columnStart[iColumn];
2288	if (model->logLevel()==63) {
2289	printf("column %d %g <= ",iColumn,columnLower[iColumn]);
2290	for (CoinBigIndex i=start;i<start+iRow;i++)
2291	printf("( %d, %g) ",row[i],element[i]);
2292	printf("<= %g\n",columnUpper[iColumn]);
2293	}
2294	}
2295	}
2296	jColumn =kColumn;
2297	}
2298	}
2299	delete [] order;
2300	delete [] other;
2301	delete [] number;
2302	}
2303	// get row copy
2304	CoinPackedMatrix rowCopy = *matrix;
2305	rowCopy.reverseOrdering();
2306	const int * rowLength = rowCopy.getVectorLengths();
2307	number = new int[numberColumns+1];
2308	memset(number,0,(numberColumns+1)*sizeof(int));
2309	for (iRow=0;iRow<numberRows;iRow++) {
2310	int length=rowLength[iRow];
2311	number[length]++;
2312	}
2313	if (number[0])
2314	printf("** %d rows have no entries\n",number[0]);
2315	k=0;
2316	for (iColumn=1;iColumn<=numberColumns;iColumn++) {
2317	if (number[iColumn]) {
2318	k++;
2319	printf("%d rows have %d entries\n",number[iColumn],iColumn);
2320	if (k==kMax)
2321	break;
2322	}
2323	}
2324	if (k<numberColumns) {
2325	int kk=k;
2326	k=0;
2327	for (iColumn=numberColumns;iColumn>=1;iColumn--) {
2328	if (number[iColumn]) {
2329	k++;
2330	if (k==kMax)
2331	break;
2332	}
2333	}
2334	if (k>kk) {
2335	printf("\n .........\n\n");
2336	iColumn=k;
2337	k=0;
2338	for (;iColumn<numberColumns;iColumn++) {
2339	if (number[iColumn]) {
2340	k++;
2341	printf("%d rows have %d entries\n",number[iColumn],iColumn);
2342	if (k==kMax)
2343	break;
2344	}
2345	}
2346	}
2347	}
2348	if (model->logLevel()==63
2349	#ifdef SYM
2350	\|\|true
2351	#endif
2352	) {
2353	int * column = rowCopy.getMutableIndices();
2354	const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
2355	double * element = rowCopy.getMutableElements();
2356	int * order = new int[numberRows];
2357	int * other = new int[numberRows];
2358	for (iRow=0;iRow<numberRows;iRow++) {
2359	int length=rowLength[iRow];
2360	order[iRow]=iRow;
2361	other[iRow]=length;
2362	CoinBigIndex start = rowStart[iRow];
2363	CoinSort_2(column+start,column+start+length,element+start);
2364	}
2365	CoinSort_2(other,other+numberRows,order);
2366	int jRow=number[0]+number[1];
2367	double * weight = new double[numberRows];
2368	double * randomColumn = new double[numberColumns+1];
2369	double * randomRow = new double [numberRows+1];
2370	int * sortRow = new int [numberRows];
2371	int * possibleRow = new int [numberRows];
2372	int * backRow = new int [numberRows];
2373	int * stackRow = new int [numberRows];
2374	int * sortColumn = new int [numberColumns];
2375	int * possibleColumn = new int [numberColumns];
2376	int * backColumn = new int [numberColumns];
2377	int * backColumn2 = new int [numberColumns];
2378	int * mapRow = new int [numberRows];
2379	int * mapColumn = new int [numberColumns];
2380	int * stackColumn = new int [numberColumns];
2381	double randomLower = CoinDrand48();
2382	double randomUpper = CoinDrand48();
2383	double randomInteger = CoinDrand48();
2384	int * startAdd = new int[numberRows+1];
2385	int * columnAdd = new int [2*numberElements];
2386	double * elementAdd = new double[2*numberElements];
2387	int nAddRows=0;
2388	startAdd[0]=0;
2389	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2390	randomColumn[iColumn] = CoinDrand48();
2391	backColumn2[iColumn]=-1;
2392	}
2393	for (iColumn=2;iColumn<=numberColumns;iColumn++) {
2394	if (number[iColumn]) {
2395	printf("XX %d rows have %d entries\n",number[iColumn],iColumn);
2396	int kRow = jRow+number[iColumn];
2397	sortOnOther(column,rowStart,
2398	order+jRow,other,number[iColumn],iColumn,0);
2399	// Now print etc
2400	if (iColumn<500000) {
2401	int nLook=0;
2402	for (int lRow =jRow;lRow<kRow;lRow++) {
2403	iRow = order[lRow];
2404	CoinBigIndex start = rowStart[iRow];
2405	if (model->logLevel()==63) {
2406	printf("row %d %g <= ",iRow,rowLower[iRow]);
2407	for (CoinBigIndex i=start;i<start+iColumn;i++)
2408	printf("( %d, %g) ",column[i],element[i]);
2409	printf("<= %g\n",rowUpper[iRow]);
2410	}
2411	int first = column[start];
2412	double sum=0.0;
2413	for (CoinBigIndex i=start;i<start+iColumn;i++) {
2414	int jColumn = column[i];
2415	double value = element[i];
2416	jColumn -= first;
2417	assert (jColumn>=0);
2418	sum += value*randomColumn[jColumn];
2419	}
2420	if (rowLower[iRow]>-1.0e30&&rowLower[iRow])
2421	sum += rowLower[iRow]*randomLower;
2422	else if (!rowLower[iRow])
2423	sum += 1.234567e-7*randomLower;
2424	if (rowUpper[iRow]<1.0e30&&rowUpper[iRow])
2425	sum += rowUpper[iRow]*randomUpper;
2426	else if (!rowUpper[iRow])
2427	sum += 1.234567e-7*randomUpper;
2428	sortRow[nLook]=iRow;
2429	randomRow[nLook++]=sum;
2430	// best way is to number unique elements and bounds and use
2431	if (fabs(sum)>1.0e4)
2432	sum *= 1.0e-6;
2433	weight[iRow]=sum;
2434	}
2435	assert (nLook<=numberRows);
2436	CoinSort_2(randomRow,randomRow+nLook,sortRow);
2437	randomRow[nLook]=COIN_DBL_MAX;
2438	double last=-COIN_DBL_MAX;
2439	int iLast=-1;
2440	for (int iLook=0;iLook<nLook+1;iLook++) {
2441	if (randomRow[iLook]>last) {
2442	if (iLast>=0) {
2443	int n=iLook-iLast;
2444	if (n>1) {
2445	//printf("%d rows possible?\n",n);
2446	}
2447	}
2448	iLast=iLook;
2449	last=randomRow[iLook];
2450	}
2451	}
2452	}
2453	jRow =kRow;
2454	}
2455	}
2456	CoinPackedMatrix columnCopy = *matrix;
2457	const int * columnLength = columnCopy.getVectorLengths();
2458	const int * row = columnCopy.getIndices();
2459	const CoinBigIndex * columnStart = columnCopy.getVectorStarts();
2460	const double * elementByColumn = columnCopy.getElements();
2461	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2462	int length=columnLength[iColumn];
2463	CoinBigIndex start = columnStart[iColumn];
2464	double sum = objective[iColumn];
2465	if (columnLower[iColumn]>-1.0e30&&columnLower[iColumn])
2466	sum += columnLower[iColumn]*randomLower;
2467	else if (!columnLower[iColumn])
2468	sum += 1.234567e-7*randomLower;
2469	if (columnUpper[iColumn]<1.0e30&&columnUpper[iColumn])
2470	sum += columnUpper[iColumn]*randomUpper;
2471	else if (!columnUpper[iColumn])
2472	sum += 1.234567e-7*randomUpper;
2473	if (model->isInteger(iColumn))
2474	sum += 9.87654321e-6*randomInteger;
2475	for (CoinBigIndex i=start;i<start+length;i++) {
2476	int iRow = row[i];
2477	sum += elementByColumn[i]*weight[iRow];
2478	}
2479	sortColumn[iColumn]=iColumn;
2480	randomColumn[iColumn]=sum;
2481	}
2482	{
2483	CoinSort_2(randomColumn,randomColumn+numberColumns,sortColumn);
2484	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2485	int i=sortColumn[iColumn];
2486	backColumn[i]=iColumn;
2487	}
2488	randomColumn[numberColumns]=COIN_DBL_MAX;
2489	double last=-COIN_DBL_MAX;
2490	int iLast=-1;
2491	for (int iLook=0;iLook<numberColumns+1;iLook++) {
2492	if (randomColumn[iLook]>last) {
2493	if (iLast>=0) {
2494	int n=iLook-iLast;
2495	if (n>1) {
2496	//printf("%d columns possible?\n",n);
2497	}
2498	for (int i=iLast;i<iLook;i++) {
2499	possibleColumn[sortColumn[i]]=n;
2500	}
2501	}
2502	iLast=iLook;
2503	last=randomColumn[iLook];
2504	}
2505	}
2506	for (iRow =0;iRow<numberRows;iRow++) {
2507	CoinBigIndex start = rowStart[iRow];
2508	double sum=0.0;
2509	int length=rowLength[iRow];
2510	for (CoinBigIndex i=start;i<start+length;i++) {
2511	int jColumn = column[i];
2512	double value = element[i];
2513	jColumn=backColumn[jColumn];
2514	sum += value*randomColumn[jColumn];
2515	//if (iColumn==23089\|\|iRow==23729)
2516	//printf("row %d cola %d colb %d value %g rand %g sum %g\n",
2517	// iRow,jColumn,column[i],value,randomColumn[jColumn],sum);
2518	}
2519	sortRow[iRow]=iRow;
2520	randomRow[iRow]=weight[iRow];
2521	randomRow[iRow]=sum;
2522	}
2523	CoinSort_2(randomRow,randomRow+numberRows,sortRow);
2524	for (iRow=0;iRow<numberRows;iRow++) {
2525	int i=sortRow[iRow];
2526	backRow[i]=iRow;
2527	}
2528	randomRow[numberRows]=COIN_DBL_MAX;
2529	last=-COIN_DBL_MAX;
2530	iLast=-1;
2531	// Do backward indices from order
2532	for (iRow=0;iRow<numberRows;iRow++) {
2533	other[order[iRow]]=iRow;
2534	}
2535	for (int iLook=0;iLook<numberRows+1;iLook++) {
2536	if (randomRow[iLook]>last) {
2537	if (iLast>=0) {
2538	int n=iLook-iLast;
2539	if (n>1) {
2540	//printf("%d rows possible?\n",n);
2541	// Within group sort as for original "order"
2542	for (int i=iLast;i<iLook;i++) {
2543	int jRow=sortRow[i];
2544	order[i]=other[jRow];
2545	}
2546	CoinSort_2(order+iLast,order+iLook,sortRow+iLast);
2547	}
2548	for (int i=iLast;i<iLook;i++) {
2549	possibleRow[sortRow[i]]=n;
2550	}
2551	}
2552	iLast=iLook;
2553	last=randomRow[iLook];
2554	}
2555	}
2556	// Temp out
2557	for (int iLook=0;iLook<numberRows-1000000;iLook++) {
2558	iRow=sortRow[iLook];
2559	CoinBigIndex start = rowStart[iRow];
2560	int length=rowLength[iRow];
2561	int numberPossible = possibleRow[iRow];
2562	for (CoinBigIndex i=start;i<start+length;i++) {
2563	int jColumn = column[i];
2564	if (possibleColumn[jColumn]!=numberPossible)
2565	numberPossible=-1;
2566	}
2567	int n=numberPossible;
2568	if (numberPossible>1) {
2569	//printf("pppppossible %d\n",numberPossible);
2570	for (int jLook=iLook+1;jLook<iLook+numberPossible;jLook++) {
2571	int jRow=sortRow[jLook];
2572	CoinBigIndex start2 = rowStart[jRow];
2573	assert (numberPossible==possibleRow[jRow]);
2574	assert(length==rowLength[jRow]);
2575	for (CoinBigIndex i=start2;i<start2+length;i++) {
2576	int jColumn = column[i];
2577	if (possibleColumn[jColumn]!=numberPossible)
2578	numberPossible=-1;
2579	}
2580	}
2581	if (numberPossible<2) {
2582	// switch off
2583	for (int jLook=iLook;jLook<iLook+n;jLook++)
2584	possibleRow[sortRow[jLook]]=-1;
2585	}
2586	// skip rest
2587	iLook += n-1;
2588	} else {
2589	possibleRow[iRow]=-1;
2590	}
2591	}
2592	for (int iLook=0;iLook<numberRows;iLook++) {
2593	iRow=sortRow[iLook];
2594	int numberPossible = possibleRow[iRow];
2595	// Only if any integers
2596	int numberIntegers=0;
2597	CoinBigIndex start = rowStart[iRow];
2598	int length=rowLength[iRow];
2599	for (CoinBigIndex i=start;i<start+length;i++) {
2600	int jColumn = column[i];
2601	if (model->isInteger(jColumn))
2602	numberIntegers++;
2603	}
2604	if (numberPossible>1&&!numberIntegers) {
2605	//printf("possible %d - but no integers\n",numberPossible);
2606	}
2607	if (numberPossible>1&&(numberIntegers\|\|true)) {
2608	//
2609	printf("possible %d - %d integers\n",numberPossible,numberIntegers);
2610	int lastLook=iLook;
2611	int nMapRow=-1;
2612	for (int jLook=iLook+1;jLook<iLook+numberPossible;jLook++) {
2613	// stop if too many failures
2614	if (jLook>iLook+10&&nMapRow<0)
2615	break;
2616	// Create identity mapping
2617	int i;
2618	for (i=0;i<numberRows;i++)
2619	mapRow[i]=i;
2620	for (i=0;i<numberColumns;i++)
2621	mapColumn[i]=i;
2622	int offset=jLook-iLook;
2623	int nStackC=0;
2624	// build up row and column mapping
2625	int nStackR=1;
2626	stackRow[0]=iLook;
2627	bool good=true;
2628	while (nStackR) {
2629	nStackR--;
2630	int look1 = stackRow[nStackR];
2631	int look2 = look1 + offset;
2632	assert (randomRow[look1]==randomRow[look2]);
2633	int row1=sortRow[look1];
2634	int row2=sortRow[look2];
2635	assert (mapRow[row1]==row1);
2636	assert (mapRow[row2]==row2);
2637	mapRow[row1]=row2;
2638	mapRow[row2]=row1;
2639	CoinBigIndex start1 = rowStart[row1];
2640	CoinBigIndex offset2 = rowStart[row2]-start1;
2641	int length=rowLength[row1];
2642	assert( length==rowLength[row2]);
2643	for (CoinBigIndex i=start1;i<start1+length;i++) {
2644	int jColumn1 = column[i];
2645	int jColumn2 = column[i+offset2];
2646	if (randomColumn[backColumn[jColumn1]]!=
2647	randomColumn[backColumn[jColumn2]]) {
2648	good=false;
2649	break;
2650	}
2651	if (mapColumn[jColumn1]==jColumn1) {
2652	// not touched
2653	assert (mapColumn[jColumn2]==jColumn2);
2654	if (jColumn1!=jColumn2) {
2655	// Put on stack
2656	mapColumn[jColumn1]=jColumn2;
2657	mapColumn[jColumn2]=jColumn1;
2658	stackColumn[nStackC++]=jColumn1;
2659	}
2660	} else {
2661	if(mapColumn[jColumn1]!=jColumn2\|\|
2662	mapColumn[jColumn2]!=jColumn1) {
2663	// bad
2664	good=false;
2665	printf("bad col\n");
2666	break;
2667	}
2668	}
2669	}
2670	if (!good)
2671	break;
2672	while (nStackC) {
2673	nStackC--;
2674	int iColumn = stackColumn[nStackC];
2675	int iColumn2 = mapColumn[iColumn];
2676	assert (iColumn!=iColumn2);
2677	int length=columnLength[iColumn];
2678	assert (length==columnLength[iColumn2]);
2679	CoinBigIndex start = columnStart[iColumn];
2680	CoinBigIndex offset2 = columnStart[iColumn2]-start;
2681	for (CoinBigIndex i=start;i<start+length;i++) {
2682	int iRow = row[i];
2683	int iRow2 = row[i+offset2];
2684	if (mapRow[iRow]==iRow) {
2685	// First (but be careful)
2686	if (iRow!=iRow2) {
2687	//mapRow[iRow]=iRow2;
2688	//mapRow[iRow2]=iRow;
2689	int iBack=backRow[iRow];
2690	int iBack2=backRow[iRow2];
2691	if (randomRow[iBack]==randomRow[iBack2]&&
2692	iBack2-iBack==offset) {
2693	stackRow[nStackR++]=iBack;
2694	} else {
2695	//printf("randomRow diff - weights %g %g\n",
2696	// weight[iRow],weight[iRow2]);
2697	// bad
2698	good=false;
2699	break;
2700	}
2701	}
2702	} else {
2703	if(mapRow[iRow]!=iRow2\|\|
2704	mapRow[iRow2]!=iRow) {
2705	// bad
2706	good=false;
2707	printf("bad row\n");
2708	break;
2709	}
2710	}
2711	}
2712	if (!good)
2713	break;
2714	}
2715	}
2716	// then check OK
2717	if (good) {
2718	for (iRow =0;iRow<numberRows;iRow++) {
2719	CoinBigIndex start = rowStart[iRow];
2720	int length=rowLength[iRow];
2721	if (mapRow[iRow]==iRow) {
2722	for (CoinBigIndex i=start;i<start+length;i++) {
2723	int jColumn = column[i];
2724	backColumn2[jColumn]=i-start;
2725	}
2726	for (CoinBigIndex i=start;i<start+length;i++) {
2727	int jColumn = column[i];
2728	if (mapColumn[jColumn]!=jColumn) {
2729	int jColumn2 =mapColumn[jColumn];
2730	CoinBigIndex i2 = backColumn2[jColumn2];
2731	if (i2<0) {
2732	good=false;
2733	} else if (element[i]!=element[i2+start]) {
2734	good=false;
2735	}
2736	}
2737	}
2738	for (CoinBigIndex i=start;i<start+length;i++) {
2739	int jColumn = column[i];
2740	backColumn2[jColumn]=-1;
2741	}
2742	} else {
2743	int row2 = mapRow[iRow];
2744	assert (iRow = mapRow[row2]);
2745	if (rowLower[iRow]!=rowLower[row2]\|\|
2746	rowLower[row2]!=rowLower[iRow])
2747	good=false;
2748	CoinBigIndex offset2 = rowStart[row2]-start;
2749	for (CoinBigIndex i=start;i<start+length;i++) {
2750	int jColumn = column[i];
2751	double value = element[i];
2752	int jColumn2 = column[i+offset2];
2753	double value2 = element[i+offset2];
2754	if (value!=value2\|\|mapColumn[jColumn]!=jColumn2\|\|
2755	mapColumn[jColumn2]!=jColumn)
2756	good=false;
2757	}
2758	}
2759	}
2760	if (good) {
2761	// check rim
2762	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2763	if (mapColumn[iColumn]!=iColumn) {
2764	int iColumn2 = mapColumn[iColumn];
2765	if(objective[iColumn]!=objective[iColumn2])
2766	good=false;
2767	if (columnLower[iColumn]!=columnLower[iColumn2])
2768	good=false;
2769	if (columnUpper[iColumn]!=columnUpper[iColumn2])
2770	good=false;
2771	if (model->isInteger(iColumn)!=model->isInteger(iColumn2))
2772	good=false;
2773	}
2774	}
2775	}
2776	if (good) {
2777	// temp
2778	if (nMapRow<0) {
2779	//const double * solution = model->primalColumnSolution();
2780	// find mapped
2781	int nMapColumn=0;
2782	for (int i=0;i<numberColumns;i++) {
2783	if (mapColumn[i]>i)
2784	nMapColumn++;
2785	}
2786	nMapRow=0;
2787	int kRow=-1;
2788	for (int i=0;i<numberRows;i++) {
2789	if (mapRow[i]>i) {
2790	nMapRow++;
2791	kRow=i;
2792	}
2793	}
2794	printf("%d columns, %d rows\n",nMapColumn,nMapRow);
2795	if (nMapRow==1) {
2796	CoinBigIndex start = rowStart[kRow];
2797	int length=rowLength[kRow];
2798	printf("%g <= ",rowLower[kRow]);
2799	for (CoinBigIndex i=start;i<start+length;i++) {
2800	int jColumn = column[i];
2801	if (mapColumn[jColumn]!=jColumn)
2802	printf("* ");
2803	printf("%d,%g ",jColumn,element[i]);
2804	}
2805	printf("<= %g\n",rowUpper[kRow]);
2806	}
2807	}
2808	// temp
2809	int row1=sortRow[lastLook];
2810	int row2=sortRow[jLook];
2811	lastLook=jLook;
2812	CoinBigIndex start1 = rowStart[row1];
2813	CoinBigIndex offset2 = rowStart[row2]-start1;
2814	int length=rowLength[row1];
2815	assert( length==rowLength[row2]);
2816	CoinBigIndex put=startAdd[nAddRows];
2817	double multiplier=length<11 ? 2.0 : 1.125;
2818	double value=1.0;
2819	for (CoinBigIndex i=start1;i<start1+length;i++) {
2820	int jColumn1 = column[i];
2821	int jColumn2 = column[i+offset2];
2822	columnAdd[put]=jColumn1;
2823	elementAdd[put++]=value;
2824	columnAdd[put]=jColumn2;
2825	elementAdd[put++]=-value;
2826	value *= multiplier;
2827	}
2828	nAddRows++;
2829	startAdd[nAddRows]=put;
2830	} else {
2831	printf("ouch - did not check out as good\n");
2832	}
2833	}
2834	}
2835	// skip rest
2836	iLook += numberPossible-1;
2837	}
2838	}
2839	}
2840	if (nAddRows) {
2841	double * lower = new double [nAddRows];
2842	double * upper = new double[nAddRows];
2843	int i;
2844	//const double * solution = model->primalColumnSolution();
2845	for (i=0;i<nAddRows;i++) {
2846	lower[i]=0.0;
2847	upper[i]=COIN_DBL_MAX;
2848	}
2849	printf("Adding %d rows with %d elements\n",nAddRows,
2850	startAdd[nAddRows]);
2851	//ClpSimplex newModel(*model);
2852	//newModel.addRows(nAddRows,lower,upper,startAdd,columnAdd,elementAdd);
2853	//newModel.writeMps("modified.mps");
2854	delete [] lower;
2855	delete [] upper;
2856	}
2857	delete [] startAdd;
2858	delete [] columnAdd;
2859	delete [] elementAdd;
2860	delete [] order;
2861	delete [] other;
2862	delete [] randomColumn;
2863	delete [] weight;
2864	delete [] randomRow;
2865	delete [] sortRow;
2866	delete [] backRow;
2867	delete [] possibleRow;
2868	delete [] sortColumn;
2869	delete [] backColumn;
2870	delete [] backColumn2;
2871	delete [] possibleColumn;
2872	delete [] mapRow;
2873	delete [] mapColumn;
2874	delete [] stackRow;
2875	delete [] stackColumn;
2876	}
2877	delete [] number;
2878	// Now do breakdown of ranges
2879	breakdown("Elements",numberElements,elementByColumn);
2880	breakdown("RowLower",numberRows,rowLower);
2881	breakdown("RowUpper",numberRows,rowUpper);
2882	breakdown("ColumnLower",numberColumns,columnLower);
2883	breakdown("ColumnUpper",numberColumns,columnUpper);
2884	breakdown("Objective",numberColumns,objective);
2885	}
2886	static bool maskMatches(const int * starts, char ** masks,
2887	std::string & check)
2888	{
2889	// back to char as I am old fashioned
2890	const char * checkC = check.c_str();
2891	int length = strlen(checkC);
2892	while (checkC[length-1]==' ')
2893	length--;
2894	for (int i=starts[length];i<starts[length+1];i++) {
2895	char * thisMask = masks[i];
2896	int k;
2897	for ( k=0;k<length;k++) {
2898	if (thisMask[k]!='?'&&thisMask[k]!=checkC[k])
2899	break;
2900	}
2901	if (k==length)
2902	return true;
2903	}
2904	return false;
2905	}
2906	static void clean(char * temp)
2907	{
2908	char * put = temp;
2909	while (*put>=' ')
2910	put++;
2911	*put='\0';
2912	}
2913	static void generateCode(const char * fileName,int type)
2914	{
2915	FILE * fp = fopen(fileName,"r");
2916	assert (fp);
2917	int numberLines=0;
2918	#define MAXLINES 500
2919	#define MAXONELINE 200
2920	char line[MAXLINES][MAXONELINE];
2921	while (fgets(line[numberLines],MAXONELINE,fp)) {
2922	assert (numberLines<MAXLINES);
2923	clean(line[numberLines]);
2924	numberLines++;
2925	}
2926	fclose(fp);
2927	// add in actual solve
2928	strcpy(line[numberLines],"5 clpModel->initialSolve(clpSolve);");
2929	numberLines++;
2930	fp = fopen(fileName,"w");
2931	assert (fp);
2932	char apo='"';
2933	char backslash = '\\';
2934
2935	fprintf(fp,"#include %cClpSimplex.hpp%c\n",apo,apo);
2936	fprintf(fp,"#include %cClpSolve.hpp%c\n",apo,apo);
2937
2938	fprintf(fp,"\nint main (int argc, const char *argv[])\n{\n");
2939	fprintf(fp," ClpSimplex model;\n");
2940	fprintf(fp," int status=1;\n");
2941	fprintf(fp," if (argc<2)\n");
2942	fprintf(fp," fprintf(stderr,%cPlease give file name%cn%c);\n",
2943	apo,backslash,apo);
2944	fprintf(fp," else\n");
2945	fprintf(fp," status=model.readMps(argv[1],true);\n");
2946	fprintf(fp," if (status) {\n");
2947	fprintf(fp," fprintf(stderr,%cBad readMps %%s%cn%c,argv[1]);\n",
2948	apo,backslash,apo);
2949	fprintf(fp," exit(1);\n");
2950	fprintf(fp," }\n\n");
2951	fprintf(fp," // Now do requested saves and modifications\n");
2952	fprintf(fp," ClpSimplex * clpModel = & model;\n");
2953	int wanted[9];
2954	memset(wanted,0,sizeof(wanted));
2955	wanted[0]=wanted[3]=wanted[5]=wanted[8]=1;
2956	if (type>0)
2957	wanted[1]=wanted[6]=1;
2958	if (type>1)
2959	wanted[2]=wanted[4]=wanted[7]=1;
2960	std::string header[9]=
2961	{ "","Save values","Redundant save of default values","Set changed values",
2962	"Redundant set default values","Solve","Restore values","Redundant restore values","Add to model"};
2963	for (int iType=0;iType<9;iType++) {
2964	if (!wanted[iType])
2965	continue;
2966	int n=0;
2967	int iLine;
2968	for (iLine=0;iLine<numberLines;iLine++) {
2969	if (line[iLine][0]=='0'+iType) {
2970	if (!n)
2971	fprintf(fp,"\n // %s\n\n",header[iType].c_str());
2972	n++;
2973	fprintf(fp,"%s\n",line[iLine]+1);
2974	}
2975	}
2976	}
2977	fprintf(fp,"\n // Now you would use solution etc etc\n\n");
2978	fprintf(fp," return 0;\n}\n");
2979	fclose(fp);
2980	printf("C++ file written to %s\n",fileName);
2981	}
2982	/*
2983	Version 1.00.00 October 13 2004.
2984	1.00.01 October 18. Added basis handling helped/prodded by Thorsten Koch.
2985	Also modifications to make faster with sbb (I hope I haven't broken anything).
2986	1.00.02 March 21 2005. Redid ClpNonLinearCost to save memory also redid
2987	createRim to try and improve cache characteristics.
2988	1.00.03 April 8 2005. Added Volume algorithm as crash and made code more
2989	robust on testing. Also added "either" and "tune" algorithm.
2990	1.01.01 April 12 2005. Decided to go to different numbering. Backups will
2991	be last 2 digits while middle 2 are for improvements. Still take a long
2992	time to get to 2.00.01
2993	1.01.02 May 4 2005. Will be putting in many changes - so saving stable version
2994	1.02.01 May 6 2005. Lots of changes to try and make faster and more stable in
2995	branch and cut.
2996	1.02.02 May 19 2005. Stuff for strong branching and some improvements to simplex
2997	1.03.01 May 24 2006. Lots done but I can't remember what!
2998	1.03.03 June 13 2006. For clean up after dual perturbation
2999	1.04.01 June 26 2007. Lots of changes but I got lazy
3000	1.05.00 June 27 2007. This is trunk so when gets to stable will be 1.5
3001	*/

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

Download in other formats: