Context Navigation

ClpSolve.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 3 months ago
formatting
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 280.1 KB

Line
1	/* $Id: ClpSolve.cpp 2385 2019-01-06 19:43:06Z unxusr $ */
2	// Copyright (C) 2003, 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	// This file has higher level solve functions
7
8	#include "CoinPragma.hpp"
9	#include "ClpConfig.h"
10
11	// check already here if COIN_HAS_GLPK is defined, since we do not want to get confused by a COIN_HAS_GLPK in config_coinutils.h
12	#if defined(COIN_HAS_AMD) \|\| defined(COIN_HAS_CHOLMOD) \|\| defined(COIN_HAS_GLPK)
13	#define UFL_BARRIER
14	#endif
15
16	#include <math.h>
17	#ifdef _MSC_VER
18	#include <windows.h> // for Sleep()
19	#ifdef small
20	#undef small
21	#endif
22	#else
23	#include <unistd.h> // for usleep()
24	#endif
25
26	#include "CoinHelperFunctions.hpp"
27	#include "ClpHelperFunctions.hpp"
28	#include "CoinSort.hpp"
29	#include "ClpFactorization.hpp"
30	#include "ClpSimplex.hpp"
31	#include "ClpSimplexOther.hpp"
32	#include "ClpSimplexDual.hpp"
33	#ifndef SLIM_CLP
34	#include "ClpQuadraticObjective.hpp"
35	#include "ClpInterior.hpp"
36	#include "ClpCholeskyDense.hpp"
37	#include "ClpCholeskyBase.hpp"
38	#include "ClpPlusMinusOneMatrix.hpp"
39	#include "ClpNetworkMatrix.hpp"
40	#endif
41	#include "ClpEventHandler.hpp"
42	#include "ClpLinearObjective.hpp"
43	#include "ClpSolve.hpp"
44	#include "ClpPackedMatrix.hpp"
45	#include "ClpMessage.hpp"
46	#include "CoinTime.hpp"
47	#if CLP_HAS_ABC
48	#include "CoinAbcCommon.hpp"
49	#endif
50	#ifdef ABC_INHERIT
51	#include "AbcSimplex.hpp"
52	#include "AbcSimplexFactorization.hpp"
53	#endif
54	#include "CoinStructuredModel.hpp"
55	double zz_slack_value = 0.0;
56	#ifdef CLP_USEFUL_PRINTOUT
57	double debugDouble[10];
58	int debugInt[24];
59	#endif
60
61	#include "ClpPresolve.hpp"
62	#ifndef SLIM_CLP
63	#include "Idiot.hpp"
64	#ifdef COIN_HAS_WSMP
65	#include "ClpCholeskyWssmp.hpp"
66	#include "ClpCholeskyWssmpKKT.hpp"
67	#endif
68	#include "ClpCholeskyUfl.hpp"
69	#ifdef TAUCS_BARRIER
70	#include "ClpCholeskyTaucs.hpp"
71	#endif
72	#ifdef PARDISO_BARRIER
73	#include "ClpCholeskyPardiso.hpp"
74	#endif
75	#include "ClpCholeskyMumps.hpp"
76	#ifdef COIN_HAS_VOL
77	#include "VolVolume.hpp"
78	#include "CoinHelperFunctions.hpp"
79	#include "CoinPackedMatrix.hpp"
80	#include "CoinMpsIO.hpp"
81	//#############################################################################
82
83	class lpHook : public VOL_user_hooks {
84	private:
85	lpHook(const lpHook &);
86	lpHook &operator=(const lpHook &);
87
88	private:
89	/// Pointer to dense vector of structural variable upper bounds
90	double *colupper_;
91	/// Pointer to dense vector of structural variable lower bounds
92	double *collower_;
93	/// Pointer to dense vector of objective coefficients
94	double *objcoeffs_;
95	/// Pointer to dense vector of right hand sides
96	double *rhs_;
97	/// Pointer to dense vector of senses
98	char *sense_;
99
100	/// The problem matrix in a row ordered form
101	CoinPackedMatrix rowMatrix_;
102	/// The problem matrix in a column ordered form
103	CoinPackedMatrix colMatrix_;
104
105	public:
106	lpHook(double clb, double cub, double *obj,
107	double rhs, char sense, const CoinPackedMatrix &mat);
108	virtual ~lpHook();
109
110	public:
111	// for all hooks: return value of -1 means that volume should quit
112	/** compute reduced costs
113	@param u (IN) the dual variables
114	@param rc (OUT) the reduced cost with respect to the dual values
115	*/
116	virtual int compute_rc(const VOL_dvector &u, VOL_dvector &rc);
117
118	/** Solve the subproblem for the subgradient step.
119	@param dual (IN) the dual variables
120	@param rc (IN) the reduced cost with respect to the dual values
121	@param lcost (OUT) the lagrangean cost with respect to the dual values
122	@param x (OUT) the primal result of solving the subproblem
123	@param v (OUT) b-Ax for the relaxed constraints
124	@param pcost (OUT) the primal objective value of <code>x</code>
125	*/
126	virtual int solve_subproblem(const VOL_dvector &dual, const VOL_dvector &rc,
127	double &lcost, VOL_dvector &x, VOL_dvector &v,
128	double &pcost);
129	/** Starting from the primal vector x, run a heuristic to produce
130	an integer solution
131	@param x (IN) the primal vector
132	@param heur_val (OUT) the value of the integer solution (return
133	<code>DBL_MAX</code> here if no feas sol was found
134	*/
135	virtual int heuristics(const VOL_problem &p,
136	const VOL_dvector &x, double &heur_val)
137	{
138	return 0;
139	}
140	};
141
142	//#############################################################################
143
144	lpHook::lpHook(double clb, double cub, double *obj,
145	double rhs, char sense,
146	const CoinPackedMatrix &mat)
147	{
148	colupper_ = cub;
149	collower_ = clb;
150	objcoeffs_ = obj;
151	rhs_ = rhs;
152	sense_ = sense;
153	assert(mat.isColOrdered());
154	colMatrix_.copyOf(mat);
155	rowMatrix_.reverseOrderedCopyOf(mat);
156	}
157
158	//-----------------------------------------------------------------------------
159
160	lpHook::~lpHook()
161	{
162	}
163
164	//#############################################################################
165
166	int lpHook::compute_rc(const VOL_dvector &u, VOL_dvector &rc)
167	{
168	rowMatrix_.transposeTimes(u.v, rc.v);
169	const int psize = rowMatrix_.getNumCols();
170
171	for (int i = 0; i < psize; ++i)
172	rc[i] = objcoeffs_[i] - rc[i];
173	return 0;
174	}
175
176	//-----------------------------------------------------------------------------
177
178	int lpHook::solve_subproblem(const VOL_dvector &dual, const VOL_dvector &rc,
179	double &lcost, VOL_dvector &x, VOL_dvector &v,
180	double &pcost)
181	{
182	int i;
183	const int psize = x.size();
184	const int dsize = v.size();
185
186	// compute the lagrangean solution corresponding to the reduced costs
187	for (i = 0; i < psize; ++i)
188	x[i] = (rc[i] >= 0.0) ? collower_[i] : colupper_[i];
189
190	// compute the lagrangean value (rhsdual + primalrc)
191	lcost = 0;
192	for (i = 0; i < dsize; ++i)
193	lcost += rhs_[i] * dual[i];
194	for (i = 0; i < psize; ++i)
195	lcost += x[i] * rc[i];
196
197	// compute the rhs - lhs
198	colMatrix_.times(x.v, v.v);
199	for (i = 0; i < dsize; ++i)
200	v[i] = rhs_[i] - v[i];
201
202	// compute the lagrangean primal objective
203	pcost = 0;
204	for (i = 0; i < psize; ++i)
205	pcost += x[i] * objcoeffs_[i];
206
207	return 0;
208	}
209
210	//#############################################################################
211	/** A quick inlined function to convert from lb/ub style constraint
212	definition to sense/rhs/range style */
213	inline void
214	convertBoundToSense(const double lower, const double upper,
215	char &sense, double &right,
216	double &range)
217	{
218	range = 0.0;
219	if (lower > -1.0e20) {
220	if (upper < 1.0e20) {
221	right = upper;
222	if (upper == lower) {
223	sense = 'E';
224	} else {
225	sense = 'R';
226	range = upper - lower;
227	}
228	} else {
229	sense = 'G';
230	right = lower;
231	}
232	} else {
233	if (upper < 1.0e20) {
234	sense = 'L';
235	right = upper;
236	} else {
237	sense = 'N';
238	right = 0.0;
239	}
240	}
241	}
242
243	static int
244	solveWithVolume(ClpSimplex *model, int numberPasses, int doIdiot)
245	{
246	VOL_problem volprob;
247	volprob.parm.gap_rel_precision = 0.00001;
248	volprob.parm.maxsgriters = 3000;
249	if (numberPasses > 3000) {
250	volprob.parm.maxsgriters = numberPasses;
251	volprob.parm.primal_abs_precision = 0.0;
252	volprob.parm.minimum_rel_ascent = 0.00001;
253	} else if (doIdiot > 0) {
254	volprob.parm.maxsgriters = doIdiot;
255	}
256	if (model->logLevel() < 2)
257	volprob.parm.printflag = 0;
258	else
259	volprob.parm.printflag = 3;
260	const CoinPackedMatrix *mat = model->matrix();
261	int psize = model->numberColumns();
262	int dsize = model->numberRows();
263	char *sense = new char[dsize];
264	double *rhs = new double[dsize];
265
266	// Set the lb/ub on the duals
267	volprob.dsize = dsize;
268	volprob.psize = psize;
269	volprob.dual_lb.allocate(dsize);
270	volprob.dual_ub.allocate(dsize);
271	int i;
272	const double *rowLower = model->rowLower();
273	const double *rowUpper = model->rowUpper();
274	for (i = 0; i < dsize; ++i) {
275	double range;
276	convertBoundToSense(rowLower[i], rowUpper[i],
277	sense[i], rhs[i], range);
278	switch (sense[i]) {
279	case 'E':
280	volprob.dual_lb[i] = -1.0e31;
281	volprob.dual_ub[i] = 1.0e31;
282	break;
283	case 'L':
284	volprob.dual_lb[i] = -1.0e31;
285	volprob.dual_ub[i] = 0.0;
286	break;
287	case 'G':
288	volprob.dual_lb[i] = 0.0;
289	volprob.dual_ub[i] = 1.0e31;
290	break;
291	default:
292	printf("Volume Algorithm can't work if there is a non ELG row\n");
293	return 1;
294	}
295	}
296	// Check bounds
297	double *saveLower = model->columnLower();
298	double *saveUpper = model->columnUpper();
299	bool good = true;
300	for (i = 0; i < psize; i++) {
301	if (saveLower[i] < -1.0e20 \|\| saveUpper[i] > 1.0e20) {
302	good = false;
303	break;
304	}
305	}
306	if (!good) {
307	saveLower = CoinCopyOfArray(model->columnLower(), psize);
308	saveUpper = CoinCopyOfArray(model->columnUpper(), psize);
309	for (i = 0; i < psize; i++) {
310	if (saveLower[i] < -1.0e20)
311	saveLower[i] = -1.0e20;
312	if (saveUpper[i] > 1.0e20)
313	saveUpper[i] = 1.0e20;
314	}
315	}
316	lpHook myHook(saveLower, saveUpper,
317	model->objective(),
318	rhs, sense, *mat);
319
320	volprob.solve(myHook, false /* no warmstart */);
321
322	if (saveLower != model->columnLower()) {
323	delete[] saveLower;
324	delete[] saveUpper;
325	}
326	//------------- extract the solution ---------------------------
327
328	//printf("Best lagrangean value: %f\n", volprob.value);
329
330	double avg = 0;
331	for (i = 0; i < dsize; ++i) {
332	switch (sense[i]) {
333	case 'E':
334	avg += CoinAbs(volprob.viol[i]);
335	break;
336	case 'L':
337	if (volprob.viol[i] < 0)
338	avg += (-volprob.viol[i]);
339	break;
340	case 'G':
341	if (volprob.viol[i] > 0)
342	avg += volprob.viol[i];
343	break;
344	}
345	}
346
347	//printf("Average primal constraint violation: %f\n", avg/dsize);
348
349	// volprob.dsol contains the dual solution (dual feasible)
350	// volprob.psol contains the primal solution
351	// (NOT necessarily primal feasible)
352	CoinMemcpyN(volprob.dsol.v, dsize, model->dualRowSolution());
353	CoinMemcpyN(volprob.psol.v, psize, model->primalColumnSolution());
354	return 0;
355	}
356	#endif
357	static ClpInterior *currentModel2 = NULL;
358	#endif
359	//#############################################################################
360	// Allow for interrupts
361	// But is this threadsafe ? (so switched off by option)
362
363	#include "CoinSignal.hpp"
364	static ClpSimplex *currentModel = NULL;
365	#ifdef ABC_INHERIT
366	static AbcSimplex *currentAbcModel = NULL;
367	#endif
368
369	extern "C" {
370	static void
371	#if defined(_MSC_VER)
372	__cdecl
373	#endif // _MSC_VER
374	signal_handler(int /whichSignal/)
375	{
376	if (currentModel != NULL)
377	currentModel->setMaximumIterations(0); // stop at next iterations
378	#ifdef ABC_INHERIT
379	if (currentAbcModel != NULL)
380	currentAbcModel->setMaximumIterations(0); // stop at next iterations
381	#endif
382	#ifndef SLIM_CLP
383	if (currentModel2 != NULL)
384	currentModel2->setMaximumBarrierIterations(0); // stop at next iterations
385	#endif
386	return;
387	}
388	}
389	#if ABC_INSTRUMENT > 1
390	int abcPricing[20];
391	int abcPricingDense[20];
392	static int trueNumberRows;
393	static int numberTypes;
394	#define MAX_TYPES 25
395	#define MAX_COUNT 20
396	#define MAX_FRACTION 101
397	static char *types[MAX_TYPES];
398	static double counts[MAX_TYPES][MAX_COUNT];
399	static double countsFraction[MAX_TYPES][MAX_FRACTION];
400	static double *currentCounts;
401	static double *currentCountsFraction;
402	static int currentType;
403	static double workMultiplier[MAX_TYPES];
404	static double work[MAX_TYPES];
405	static double currentWork;
406	static double otherWork[MAX_TYPES];
407	static int timesCalled[MAX_TYPES];
408	static int timesStarted[MAX_TYPES];
409	static int fractionDivider;
410	void instrument_initialize(int numberRows)
411	{
412	trueNumberRows = numberRows;
413	numberTypes = 0;
414	memset(counts, 0, sizeof(counts));
415	currentCounts = NULL;
416	memset(countsFraction, 0, sizeof(countsFraction));
417	currentCountsFraction = NULL;
418	memset(workMultiplier, 0, sizeof(workMultiplier));
419	memset(work, 0, sizeof(work));
420	memset(otherWork, 0, sizeof(otherWork));
421	memset(timesCalled, 0, sizeof(timesCalled));
422	memset(timesStarted, 0, sizeof(timesStarted));
423	currentType = -1;
424	fractionDivider = (numberRows + MAX_FRACTION - 2) / (MAX_FRACTION - 1);
425	}
426	void instrument_start(const char *type, int numberRowsEtc)
427	{
428	if (currentType >= 0)
429	instrument_end();
430	currentType = -1;
431	currentWork = 0.0;
432	for (int i = 0; i < numberTypes; i++) {
433	if (!strcmp(types[i], type)) {
434	currentType = i;
435	break;
436	}
437	}
438	if (currentType == -1) {
439	assert(numberTypes < MAX_TYPES);
440	currentType = numberTypes;
441	types[numberTypes++] = strdup(type);
442	}
443	currentCounts = &counts[currentType][0];
444	currentCountsFraction = &countsFraction[currentType][0];
445	timesStarted[currentType]++;
446	assert(trueNumberRows);
447	workMultiplier[currentType] += static_cast< double >(numberRowsEtc) / static_cast< double >(trueNumberRows);
448	}
449	void instrument_add(int count)
450	{
451	assert(currentType >= 0);
452	currentWork += count;
453	timesCalled[currentType]++;
454	if (count < MAX_COUNT - 1)
455	currentCounts[count]++;
456	else
457	currentCounts[MAX_COUNT - 1]++;
458	assert(count / fractionDivider >= 0 && count / fractionDivider < MAX_FRACTION);
459	currentCountsFraction[count / fractionDivider]++;
460	}
461	void instrument_do(const char *type, double count)
462	{
463	int iType = -1;
464	for (int i = 0; i < numberTypes; i++) {
465	if (!strcmp(types[i], type)) {
466	iType = i;
467	break;
468	}
469	}
470	if (iType == -1) {
471	assert(numberTypes < MAX_TYPES);
472	iType = numberTypes;
473	types[numberTypes++] = strdup(type);
474	}
475	timesStarted[iType]++;
476	otherWork[iType] += count;
477	}
478	void instrument_end()
479	{
480	work[currentType] += currentWork;
481	currentType = -1;
482	}
483	void instrument_end_and_adjust(double factor)
484	{
485	work[currentType] += currentWork * factor;
486	currentType = -1;
487	}
488	void instrument_print()
489	{
490	for (int iType = 0; iType < numberTypes; iType++) {
491	currentCounts = &counts[iType][0];
492	currentCountsFraction = &countsFraction[iType][0];
493	if (!otherWork[iType]) {
494	printf("%s started %d times, used %d times, work %g (average length %.1f) multiplier %g\n",
495	types[iType], timesStarted[iType], timesCalled[iType],
496	work[iType], work[iType] / (timesCalled[iType] + 1.0e-100), workMultiplier[iType] / (timesStarted[iType] + 1.0e-100));
497	int n = 0;
498	for (int i = 0; i < MAX_COUNT - 1; i++) {
499	if (currentCounts[i]) {
500	if (n == 5) {
501	n = 0;
502	printf("\n");
503	}
504	n++;
505	printf("(%d els,%.0f times) ", i, currentCounts[i]);
506	}
507	}
508	if (currentCounts[MAX_COUNT - 1]) {
509	if (n == 5) {
510	n = 0;
511	printf("\n");
512	}
513	n++;
514	printf("(>=%d els,%.0f times) ", MAX_COUNT - 1, currentCounts[MAX_COUNT - 1]);
515	}
516	printf("\n");
517	int largestFraction;
518	int nBig = 0;
519	for (largestFraction = MAX_FRACTION - 1; largestFraction >= 10; largestFraction--) {
520	double count = currentCountsFraction[largestFraction];
521	if (count && largestFraction > 10)
522	nBig++;
523	if (nBig > 4)
524	break;
525	}
526	int chunk = (largestFraction + 5) / 10;
527	int lo = 0;
528	for (int iChunk = 0; iChunk < largestFraction; iChunk += chunk) {
529	int hi = CoinMin(lo + chunk * fractionDivider, trueNumberRows);
530	double sum = 0.0;
531	for (int i = iChunk; i < CoinMin(iChunk + chunk, MAX_FRACTION); i++)
532	sum += currentCountsFraction[i];
533	if (sum)
534	printf("(%d-%d %.0f) ", lo, hi, sum);
535	lo = hi;
536	}
537	for (int i = lo / fractionDivider; i < MAX_FRACTION; i++) {
538	if (currentCountsFraction[i])
539	printf("(%d %.0f) ", i * fractionDivider, currentCountsFraction[i]);
540	}
541	printf("\n");
542	} else {
543	printf("%s started %d times, used %d times, work %g multiplier %g other work %g\n",
544	types[iType], timesStarted[iType], timesCalled[iType],
545	work[iType], workMultiplier[iType], otherWork[iType]);
546	}
547	free(types[iType]);
548	}
549	}
550	#endif
551	#if ABC_PARALLEL == 2
552	#ifndef FAKE_CILK
553	int number_cilk_workers = 0;
554	#include <cilk/cilk_api.h>
555	#endif
556	#endif
557	#ifdef ABC_INHERIT
558	AbcSimplex *
559	ClpSimplex::dealWithAbc(int solveType, int startUp,
560	bool interrupt)
561	{
562	bool keepAbc = startUp >= 10;
563	startUp = startUp % 10;
564	AbcSimplex *abcModel2 = NULL;
565	if (!this->abcState() \|\| !numberRows_ \|\| !numberColumns_) {
566	//this->readBasis("aaa.bas");
567	if (!solveType) {
568	this->dual(0);
569	} else if (solveType == 1) {
570	int ifValuesPass = startUp ? 1 : 0;
571	if (startUp == 3)
572	ifValuesPass = 2;
573	this->primal(ifValuesPass);
574	}
575	//this->writeBasis("a.bas",true);
576	} else {
577	abcModel2 = new AbcSimplex(*this);
578	if (interrupt)
579	currentAbcModel = abcModel2;
580	//if (abcSimplex_) {
581	// move factorization stuff
582	abcModel2->factorization()->synchronize(this->factorization(), abcModel2);
583	//}
584	//abcModel2->startPermanentArrays();
585	int crashState = abcModel2->abcState() & (256 + 512 + 1024);
586	abcModel2->setAbcState(CLP_ABC_WANTED \| crashState \| (abcModel2->abcState() & 15));
587	int ifValuesPass = startUp ? 1 : 0;
588	if (startUp == 3)
589	ifValuesPass = 2;
590	// temp
591	if (fabs(this->primalTolerance() - 1.001e-6) < 0.999e-9) {
592	int type = 1;
593	double diff = this->primalTolerance() - 1.001e-6;
594	printf("Diff %g\n", diff);
595	if (fabs(diff - 1.0e-10) < 1.0e-13)
596	type = 2;
597	else if (fabs(diff - 1.0e-11) < 1.0e-13)
598	type = 3;
599	#if 0
600	ClpSimplex * thisModel = static_cast<ClpSimplexOther *> (this)->dualOfModel(1.0,1.0);
601	if (thisModel) {
602	printf("Dual of model has %d rows and %d columns\n",
603	thisModel->numberRows(), thisModel->numberColumns());
604	thisModel->setOptimizationDirection(1.0);
605	Idiot info(*thisModel);
606	info.setStrategy(512 \| info.getStrategy());
607	// Allow for scaling
608	info.setStrategy(32 \| info.getStrategy());
609	info.setStartingWeight(1.0e3);
610	info.setReduceIterations(6);
611	info.crash(50, this->messageHandler(), this->messagesPointer(),false);
612	// make sure later that primal solution in correct place
613	// and has correct sign
614	abcModel2->setupDualValuesPass(thisModel->primalColumnSolution(),
615	thisModel->dualRowSolution(),type);
616	//thisModel->dual();
617	delete thisModel;
618	}
619	#else
620	if (!solveType) {
621	this->dual(0);
622	abcModel2->setupDualValuesPass(this->dualRowSolution(),
623	this->primalColumnSolution(), type);
624	} else if (solveType == 1) {
625	ifValuesPass = 1;
626	abcModel2->setStateOfProblem(abcModel2->stateOfProblem() \| VALUES_PASS);
627	Idiot info(*abcModel2);
628	info.setStrategy(512 \| info.getStrategy());
629	// Allow for scaling
630	info.setStrategy(32 \| info.getStrategy());
631	info.setStartingWeight(1.0e3);
632	info.setReduceIterations(6);
633	info.crash(200, abcModel2->messageHandler(), abcModel2->messagesPointer(), false);
634	//memcpy(abcModel2->primalColumnSolution(),this->primalColumnSolution(),
635	// this->numberColumns()*sizeof(double));
636	}
637	#endif
638	}
639	int numberCpu = this->abcState() & 15;
640	if (numberCpu == 9) {
641	numberCpu = 1;
642	#if ABC_PARALLEL == 2
643	#ifndef FAKE_CILK
644	if (number_cilk_workers > 1)
645	numberCpu = CoinMin(2 * number_cilk_workers, 8);
646	#endif
647	#endif
648	} else if (numberCpu == 10) {
649	// maximum
650	numberCpu = 4;
651	} else if (numberCpu == 10) {
652	// decide
653	if (abcModel2->getNumElements() < 5000)
654	numberCpu = 1;
655	#if ABC_PARALLEL == 2
656	#ifndef FAKE_CILK
657	else if (number_cilk_workers > 1)
658	numberCpu = CoinMin(2 * number_cilk_workers, 8);
659	#endif
660	#endif
661	else
662	numberCpu = 1;
663	} else {
664	#if ABC_PARALLEL == 2
665	#ifndef FAKE_CILK
666	char temp[3];
667	sprintf(temp, "%d", numberCpu);
668	__cilkrts_set_param("nworkers", temp);
669	printf("setting cilk workers to %d\n", numberCpu);
670	number_cilk_workers = numberCpu;
671
672	#endif
673	#endif
674	}
675	char line[200];
676	#if ABC_PARALLEL
677	#if ABC_PARALLEL == 2
678	#ifndef FAKE_CILK
679	if (!number_cilk_workers) {
680	number_cilk_workers = __cilkrts_get_nworkers();
681	sprintf(line, "%d cilk workers", number_cilk_workers);
682	handler_->message(CLP_GENERAL, messages_)
683	<< line
684	<< CoinMessageEol;
685	}
686	#endif
687	#endif
688	abcModel2->setParallelMode(numberCpu - 1);
689	#endif
690	//if (abcState()==3\|\|abcState()==4) {
691	//abcModel2->setMoreSpecialOptions((131072*2)\|abcModel2->moreSpecialOptions());
692	//}
693	//if (processTune>0&&processTune<8)
694	//abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()\|131072*processTune);
695	#if ABC_INSTRUMENT
696	double startTimeCpu = CoinCpuTime();
697	double startTimeElapsed = CoinGetTimeOfDay();
698	#if ABC_INSTRUMENT > 1
699	memset(abcPricing, 0, sizeof(abcPricing));
700	memset(abcPricingDense, 0, sizeof(abcPricing));
701	instrument_initialize(abcModel2->numberRows());
702	#endif
703	#endif
704	if (!solveType) {
705	abcModel2->ClpSimplex::doAbcDual();
706	} else if (solveType == 1) {
707	int saveOptions = abcModel2->specialOptions();
708	//if (startUp==2)
709	//abcModel2->setSpecialOptions(8192\|saveOptions);
710	abcModel2->ClpSimplex::doAbcPrimal(ifValuesPass);
711	abcModel2->setSpecialOptions(saveOptions);
712	}
713	#if ABC_INSTRUMENT
714	if (solveType < 2) {
715	double timeCpu = CoinCpuTime() - startTimeCpu;
716	double timeElapsed = CoinGetTimeOfDay() - startTimeElapsed;
717	sprintf(line, "Cpu time for %s (%d rows, %d columns %d elements) %g elapsed %g ratio %g - %d iterations",
718	this->problemName().c_str(), this->numberRows(), this->numberColumns(),
719	this->getNumElements(),
720	timeCpu, timeElapsed, timeElapsed ? timeCpu / timeElapsed : 1.0, abcModel2->numberIterations());
721	handler_->message(CLP_GENERAL, messages_)
722	<< line
723	<< CoinMessageEol;
724	#if ABC_INSTRUMENT > 1
725	{
726	int n;
727	n = 0;
728	for (int i = 0; i < 20; i++)
729	n += abcPricing[i];
730	printf("CCSparse pricing done %d times", n);
731	int n2 = 0;
732	for (int i = 0; i < 20; i++)
733	n2 += abcPricingDense[i];
734	if (n2)
735	printf(" and dense pricing done %d times\n", n2);
736	else
737	printf("\n");
738	n = 0;
739	printf("CCS");
740	for (int i = 0; i < 19; i++) {
741	if (abcPricing[i]) {
742	if (n == 5) {
743	n = 0;
744	printf("\nCCS");
745	}
746	n++;
747	printf("(%d els,%d times) ", i, abcPricing[i]);
748	}
749	}
750	if (abcPricing[19]) {
751	if (n == 5) {
752	n = 0;
753	printf("\nCCS");
754	}
755	n++;
756	printf("(>=19 els,%d times) ", abcPricing[19]);
757	}
758	if (n2) {
759	printf("CCD");
760	for (int i = 0; i < 19; i++) {
761	if (abcPricingDense[i]) {
762	if (n == 5) {
763	n = 0;
764	printf("\nCCD");
765	}
766	n++;
767	int k1 = (numberRows_ / 16) * i;
768	;
769	int k2 = CoinMin(numberRows_, k1 + (numberRows_ / 16) - 1);
770	printf("(%d-%d els,%d times) ", k1, k2, abcPricingDense[i]);
771	}
772	}
773	}
774	printf("\n");
775	}
776	instrument_print();
777	#endif
778	}
779	#endif
780	abcModel2->moveStatusToClp(this);
781	#if 1
782	if (!problemStatus_) {
783	double offset;
784	CoinMemcpyN(this->objectiveAsObject()->gradient(this,
785	this->primalColumnSolution(), offset, true),
786	numberColumns_, this->dualColumnSolution());
787	this->clpMatrix()->transposeTimes(-1.0,
788	this->dualRowSolution(),
789	this->dualColumnSolution());
790	memset(this->primalRowSolution(), 0, numberRows_ * sizeof(double));
791	this->clpMatrix()->times(1.0,
792	this->primalColumnSolution(),
793	this->primalRowSolution());
794	this->checkSolutionInternal();
795	if (sumDualInfeasibilities_ > 100.0 * dualTolerance_) {
796	#if CBC_USEFUL_PRINTING > 0
797	printf("internal check on duals failed %d %g\n",
798	numberDualInfeasibilities_, sumDualInfeasibilities_);
799	#endif
800	} else {
801	sumDualInfeasibilities_ = 0.0;
802	numberDualInfeasibilities_ = 0;
803	}
804	if (sumPrimalInfeasibilities_ > 100.0 * primalTolerance_) {
805	#if CBC_USEFUL_PRINTING > 0
806	printf("internal check on primals failed %d %g\n",
807	numberPrimalInfeasibilities_, sumPrimalInfeasibilities_);
808	#endif
809	} else {
810	sumPrimalInfeasibilities_ = 0.0;
811	numberPrimalInfeasibilities_ = 0;
812	}
813	problemStatus_ = 0;
814	}
815	#endif
816	//ClpModel::stopPermanentArrays();
817	this->setSpecialOptions(this->specialOptions() & ~65536);
818	#if 0
819	this->writeBasis("a.bas",true);
820	for (int i=0;i<this->numberRows();i++)
821	printf("%d %g\n",i,this->dualRowSolution()[i]);
822	this->dual();
823	this->writeBasis("b.bas",true);
824	for (int i=0;i<this->numberRows();i++)
825	printf("%d %g\n",i,this->dualRowSolution()[i]);
826	#endif
827	// switch off initialSolve flag
828	moreSpecialOptions_ &= ~16384;
829	//this->setNumberIterations(abcModel2->numberIterations()+this->numberIterations());
830	if (!keepAbc) {
831	delete abcModel2;
832	abcModel2 = NULL;
833	}
834	}
835	return abcModel2;
836	}
837	#endif
838	/** General solve algorithm which can do presolve
839	special options (bits)
840	1 - do not perturb
841	2 - do not scale
842	4 - use crash (default allslack in dual, idiot in primal)
843	8 - all slack basis in primal
844	16 - switch off interrupt handling
845	32 - do not try and make plus minus one matrix
846	64 - do not use sprint even if problem looks good
847	*/
848	int ClpSimplex::initialSolve(ClpSolve &options)
849	{
850	ClpSolve::SolveType method = options.getSolveType();
851	//ClpSolve::SolveType originalMethod=method;
852	ClpSolve::PresolveType presolve = options.getPresolveType();
853	int saveMaxIterations = maximumIterations();
854	int finalStatus = -1;
855	int numberIterations = 0;
856	double time1 = CoinCpuTime();
857	double timeX = time1;
858	double time2 = 0.0;
859	ClpMatrixBase *saveMatrix = NULL;
860	ClpObjective *savedObjective = NULL;
861	int idiotOptions = 0;
862	if (options.getSpecialOption(6))
863	idiotOptions = options.getExtraInfo(6) * 32768;
864	#ifdef CLP_USEFUL_PRINTOUT
865	debugInt[0] = numberRows();
866	debugInt[1] = numberColumns();
867	debugInt[2] = matrix()->getNumElements();
868	#endif
869	if (!objective_ \|\| !matrix_) {
870	// totally empty
871	handler_->message(CLP_EMPTY_PROBLEM, messages_)
872	<< 0
873	<< 0
874	<< 0
875	<< CoinMessageEol;
876	return -1;
877	} else if (!numberRows_ \|\| !numberColumns_ \|\| !getNumElements()) {
878	presolve = ClpSolve::presolveOff;
879	}
880	if (objective_->type() >= 2 && optimizationDirection_ == 0) {
881	// pretend linear
882	savedObjective = objective_;
883	// make up objective
884	double *obj = new double[numberColumns_];
885	for (int i = 0; i < numberColumns_; i++) {
886	double l = fabs(columnLower_[i]);
887	double u = fabs(columnUpper_[i]);
888	obj[i] = 0.0;
889	if (CoinMin(l, u) < 1.0e20) {
890	if (l < u)
891	obj[i] = 1.0 + randomNumberGenerator_.randomDouble() * 1.0e-2;
892	else
893	obj[i] = -1.0 - randomNumberGenerator_.randomDouble() * 1.0e-2;
894	}
895	}
896	objective_ = new ClpLinearObjective(obj, numberColumns_);
897	delete[] obj;
898	}
899	ClpSimplex *model2 = this;
900	bool interrupt = (options.getSpecialOption(2) == 0);
901	CoinSighandler_t saveSignal = static_cast< CoinSighandler_t >(0);
902	if (interrupt) {
903	currentModel = model2;
904	// register signal handler
905	saveSignal = signal(SIGINT, signal_handler);
906	}
907	// If no status array - set up basis
908	if (!status_)
909	allSlackBasis();
910	ClpPresolve *pinfo = new ClpPresolve();
911	pinfo->setSubstitution(options.substitution());
912	int presolveOptions = options.presolveActions();
913	bool presolveToFile = (presolveOptions & 0x40000000) != 0;
914	presolveOptions &= ~0x40000000;
915	if ((presolveOptions & 0xffffff) != 0)
916	pinfo->setPresolveActions(presolveOptions);
917	// switch off singletons to slacks
918	//pinfo->setDoSingletonColumn(false); // done by bits
919	int printOptions = options.getSpecialOption(5);
920	if ((printOptions & 1) != 0)
921	pinfo->statistics();
922	double timePresolve = 0.0;
923	double timeIdiot = 0.0;
924	double timeCore = 0.0;
925	eventHandler()->event(ClpEventHandler::presolveStart);
926	int savePerturbation = perturbation_;
927	int saveScaling = scalingFlag_;
928	#ifndef SLIM_CLP
929	#ifndef NO_RTTI
930	if (dynamic_cast< ClpNetworkMatrix * >(matrix_)) {
931	// network - switch off stuff
932	presolve = ClpSolve::presolveOff;
933	}
934	#else
935	if (matrix_->type() == 11) {
936	// network - switch off stuff
937	presolve = ClpSolve::presolveOff;
938	}
939	#endif
940	#endif
941	#ifndef CLPSOLVE_ACTIONS
942	#define CLPSOLVE_ACTIONS 2
943	#endif
944	#if CLPSOLVE_ACTIONS
945	bool wasAutomatic = (method == ClpSolve::automatic);
946	#endif
947	if (presolve != ClpSolve::presolveOff) {
948	bool costedSlacks = false;
949	#if CLP_INHERIT_MODE > 1
950	int numberPasses = 20;
951	#else
952	int numberPasses = 5;
953	#endif
954	if (presolve == ClpSolve::presolveNumber) {
955	numberPasses = options.getPresolvePasses();
956	presolve = ClpSolve::presolveOn;
957	} else if (presolve == ClpSolve::presolveNumberCost) {
958	numberPasses = options.getPresolvePasses();
959	presolve = ClpSolve::presolveOn;
960	costedSlacks = true;
961	// switch on singletons to slacks
962	pinfo->setDoSingletonColumn(true);
963	// gub stuff for testing
964	//pinfo->setDoGubrow(true);
965	}
966	#ifndef CLP_NO_STD
967	if (presolveToFile) {
968	// PreSolve to file - not fully tested
969	printf("Presolving to file - presolve.save\n");
970	pinfo->presolvedModelToFile(*this, "presolve.save", dblParam_[ClpPresolveTolerance],
971	false, numberPasses);
972	model2 = this;
973	} else {
974	#endif
975	model2 = pinfo->presolvedModel(*this, dblParam_[ClpPresolveTolerance],
976	false, numberPasses, true, costedSlacks);
977	#ifndef CLP_NO_STD
978	}
979	#endif
980	time2 = CoinCpuTime();
981	timePresolve = time2 - timeX;
982	handler_->message(CLP_INTERVAL_TIMING, messages_)
983	<< "Presolve" << timePresolve << time2 - time1
984	<< CoinMessageEol;
985	timeX = time2;
986	if (!model2) {
987	handler_->message(CLP_INFEASIBLE, messages_)
988	<< CoinMessageEol;
989	model2 = this;
990	eventHandler()->event(ClpEventHandler::presolveInfeasible);
991	problemStatus_ = pinfo->presolveStatus();
992	if (options.infeasibleReturn() \|\| (moreSpecialOptions_ & 1) != 0) {
993	delete pinfo;
994	return -1;
995	}
996	presolve = ClpSolve::presolveOff;
997	} else {
998	#if 0 //def ABC_INHERIT
999	{
1000	AbcSimplex * abcModel2=new AbcSimplex(*model2);
1001	delete model2;
1002	model2=abcModel2;
1003	pinfo->setPresolvedModel(model2);
1004	}
1005	#else
1006	//ClpModel::stopPermanentArrays();
1007	//setSpecialOptions(specialOptions()&~65536);
1008	// try setting tolerances up
1009	#if CLPSOLVE_ACTIONS
1010	bool changeTolerances = wasAutomatic;
1011	#if CLPSOLVE_ACTIONS > 1
1012	changeTolerances = true;
1013	#endif
1014	if (changeTolerances && model2 != this) {
1015	#define CLP_NEW_TOLERANCE 1.0e-6
1016	if (model2->primalTolerance() == 1.0e-7 && model2->dualTolerance() == 1.0e-7) {
1017	model2->setPrimalTolerance(CLP_NEW_TOLERANCE);
1018	model2->setDualTolerance(CLP_NEW_TOLERANCE);
1019	}
1020	}
1021	#endif
1022	#endif
1023	model2->eventHandler()->setSimplex(model2);
1024	int rcode = model2->eventHandler()->event(ClpEventHandler::presolveSize);
1025	// see if too big or small
1026	if (rcode == 2) {
1027	delete model2;
1028	delete pinfo;
1029	return -2;
1030	} else if (rcode == 3) {
1031	delete model2;
1032	delete pinfo;
1033	return -3;
1034	}
1035	}
1036	model2->setMoreSpecialOptions(model2->moreSpecialOptions() & (~1024));
1037	model2->eventHandler()->setSimplex(model2);
1038	// We may be better off using original (but if dual leave because of bounds)
1039	if (presolve != ClpSolve::presolveOff && numberRows_ < 1.01 * model2->numberRows_ && numberColumns_ < 1.01 * model2->numberColumns_
1040	&& model2 != this) {
1041	if (method != ClpSolve::useDual \|\| (numberRows_ == model2->numberRows_ && numberColumns_ == model2->numberColumns_)) {
1042	delete model2;
1043	model2 = this;
1044	presolve = ClpSolve::presolveOff;
1045	}
1046	}
1047	}
1048	#ifdef CLP_USEFUL_PRINTOUT
1049	debugInt[3] = model2->numberRows();
1050	debugInt[4] = model2->numberColumns();
1051	debugInt[5] = model2->matrix()->getNumElements();
1052	// analyze
1053	{
1054	double time1 = CoinCpuTime();
1055	int numberColumns = model2->numberColumns();
1056	const double *columnLower = model2->columnLower();
1057	const double *columnUpper = model2->columnUpper();
1058	int numberRows = model2->numberRows();
1059	const double *rowLower = model2->rowLower();
1060	const double *rowUpper = model2->rowUpper();
1061	const double *objective = model2->objective();
1062	CoinPackedMatrix *matrix = model2->matrix();
1063	CoinBigIndex numberElements = matrix->getNumElements();
1064	const int *columnLength = matrix->getVectorLengths();
1065	//const CoinBigIndex * columnStart = matrix->getVectorStarts();
1066	const double *elementByColumn = matrix->getElements();
1067	const int *row = matrix->getIndices();
1068	int *rowCount = new int[numberRows];
1069	memset(rowCount, 0, numberRows * sizeof(int));
1070	int n = CoinMax(2 * numberRows, numberElements);
1071	n = CoinMax(2 * numberColumns, n);
1072	double *check = new double[n];
1073	memcpy(check, elementByColumn, numberElements * sizeof(double));
1074	for (int i = 0; i < numberElements; i++) {
1075	check[i] = fabs(check[i]);
1076	rowCount[row[i]]++;
1077	}
1078	int largestIndex = 0;
1079	for (int i = 0; i < numberColumns; i++) {
1080	largestIndex = CoinMax(largestIndex, columnLength[i]);
1081	}
1082	debugInt[12] = largestIndex;
1083	largestIndex = 0;
1084	for (int i = 0; i < numberRows; i++) {
1085	largestIndex = CoinMax(largestIndex, rowCount[i]);
1086	}
1087	n = numberElements;
1088	delete[] rowCount;
1089	debugInt[11] = largestIndex;
1090	std::sort(check, check + n);
1091	debugDouble[4] = check[0];
1092	debugDouble[5] = check[n - 1];
1093	int nAtOne = 0;
1094	int nDifferent = 0;
1095	double last = -COIN_DBL_MAX;
1096	for (int i = 0; i < n; i++) {
1097	if (fabs(last - check[i]) > 1.0e-12) {
1098	nDifferent++;
1099	last = check[i];
1100	}
1101	if (check[i] == 1.0)
1102	nAtOne++;
1103	}
1104	debugInt[10] = nDifferent;
1105	debugInt[15] = nAtOne;
1106	int nInf = 0;
1107	int nZero = 0;
1108	n = 0;
1109	for (int i = 0; i < numberRows; i++) {
1110	double value = fabs(rowLower[i]);
1111	if (!value)
1112	nZero++;
1113	else if (value != COIN_DBL_MAX)
1114	check[n++] = value;
1115	else
1116	nInf++;
1117	}
1118	for (int i = 0; i < numberRows; i++) {
1119	double value = fabs(rowUpper[i]);
1120	if (!value)
1121	nZero++;
1122	else if (value != COIN_DBL_MAX)
1123	check[n++] = value;
1124	else
1125	nInf++;
1126	}
1127	debugInt[16] = nInf;
1128	debugInt[20] = nZero;
1129	if (n) {
1130	std::sort(check, check + n);
1131	debugDouble[0] = check[0];
1132	debugDouble[1] = check[n - 1];
1133	} else {
1134	debugDouble[0] = 0.0;
1135	debugDouble[1] = 0.0;
1136	}
1137	nAtOne = 0;
1138	nDifferent = 0;
1139	last = -COIN_DBL_MAX;
1140	for (int i = 0; i < n; i++) {
1141	if (fabs(last - check[i]) > 1.0e-12) {
1142	nDifferent++;
1143	last = check[i];
1144	}
1145	if (check[i] == 1.0)
1146	nAtOne++;
1147	}
1148	debugInt[8] = nDifferent;
1149	debugInt[13] = nAtOne;
1150	nZero = 0;
1151	n = 0;
1152	for (int i = 0; i < numberColumns; i++) {
1153	double value = fabs(objective[i]);
1154	if (value)
1155	check[n++] = value;
1156	else
1157	nZero++;
1158	}
1159	debugInt[21] = nZero;
1160	if (n) {
1161	std::sort(check, check + n);
1162	debugDouble[2] = check[0];
1163	debugDouble[3] = check[n - 1];
1164	} else {
1165	debugDouble[2] = 0.0;
1166	debugDouble[3] = 0.0;
1167	}
1168	nAtOne = 0;
1169	nDifferent = 0;
1170	last = -COIN_DBL_MAX;
1171	for (int i = 0; i < n; i++) {
1172	if (fabs(last - check[i]) > 1.0e-12) {
1173	nDifferent++;
1174	last = check[i];
1175	}
1176	if (check[i] == 1.0)
1177	nAtOne++;
1178	}
1179	debugInt[9] = nDifferent;
1180	debugInt[14] = nAtOne;
1181	nInf = 0;
1182	nZero = 0;
1183	n = 0;
1184	for (int i = 0; i < numberColumns; i++) {
1185	double value = fabs(columnLower[i]);
1186	if (!value)
1187	nZero++;
1188	else if (value != COIN_DBL_MAX)
1189	check[n++] = value;
1190	else
1191	nInf++;
1192	}
1193	for (int i = 0; i < numberColumns; i++) {
1194	double value = fabs(columnUpper[i]);
1195	if (!value)
1196	nZero++;
1197	else if (value != COIN_DBL_MAX)
1198	check[n++] = value;
1199	else
1200	nInf++;
1201	}
1202	debugInt[17] = nInf;
1203	double smallestColBound;
1204	double largestColBound;
1205	if (n) {
1206	std::sort(check, check + n);
1207	smallestColBound = check[0];
1208	largestColBound = check[n - 1];
1209	} else {
1210	smallestColBound = 0.0;
1211	largestColBound = 0.0;
1212	}
1213	nAtOne = 0;
1214	nDifferent = 0;
1215	last = -COIN_DBL_MAX;
1216	for (int i = 0; i < n; i++) {
1217	if (fabs(last - check[i]) > 1.0e-12) {
1218	nDifferent++;
1219	last = check[i];
1220	}
1221	if (check[i] == 1.0)
1222	nAtOne++;
1223	}
1224	//debugInt[8]=nDifferent;
1225	//debugInt[13]=nAtOne;
1226	printf("BENCHMARK_STATS rhs %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
1227	debugInt[8], debugDouble[0], debugDouble[1], debugInt[13], debugInt[16], debugInt[20]);
1228	printf("BENCHMARK_STATS col %d different - %g -> %g (%d at one, %d infinite, %d zero)\n",
1229	nDifferent, smallestColBound, largestColBound, nAtOne, nInf, nZero);
1230	printf("BENCHMARK_STATS els %d different - %g -> %g (%d at one) - longest r,c %d,%d\n",
1231	debugInt[10], debugDouble[4], debugDouble[5], debugInt[15],
1232	debugInt[11], debugInt[12]);
1233	printf("BENCHMARK_STATS obj %d different - %g -> %g (%d at one, %d zero) - time %g\n",
1234	debugInt[9], debugDouble[2], debugDouble[3], debugInt[14], debugInt[21],
1235	CoinCpuTime() - time1);
1236	delete[] check;
1237	}
1238	#endif
1239	if (interrupt)
1240	currentModel = model2;
1241	int saveMoreOptions = moreSpecialOptions_;
1242	// For below >0 overrides
1243	// 0 means no, -1 means maybe
1244	int doIdiot = 0;
1245	int doCrash = 0;
1246	int doSprint = 0;
1247	int doSlp = 0;
1248	int primalStartup = 1;
1249	model2->eventHandler()->event(ClpEventHandler::presolveBeforeSolve);
1250	#if CLP_POOL_MATRIX
1251	if (vectorMode() >= 10) {
1252	ClpPoolMatrix poolMatrix = new ClpPoolMatrix(model2->matrix());
1253	char output[80];
1254	int numberDifferent = poolMatrix->getNumDifferentElements();
1255	if (numberDifferent > 0) {
1256	sprintf(output, "Pool matrix has %d different values",
1257	numberDifferent);
1258	model2->replaceMatrix(poolMatrix, true);
1259	} else {
1260	delete poolMatrix;
1261	sprintf(output, "Pool matrix has more than %d different values - no good",
1262	-numberDifferent);
1263	}
1264	handler_->message(CLP_GENERAL, messages_) << output
1265	<< CoinMessageEol;
1266	}
1267	#endif
1268	int tryItSave = 0;
1269	#if CLPSOLVE_ACTIONS
1270	if (method == ClpSolve::automatic)
1271	model2->moreSpecialOptions_ \|= 8192; // stop switch over
1272	#endif
1273	// switch to primal from automatic if just one cost entry
1274	if (method == ClpSolve::automatic && model2->numberColumns() > 5000
1275	#ifndef CLPSOLVE_ACTIONS
1276	&& (specialOptions_ & 1024) != 0
1277	#endif
1278	) {
1279	// look at original model for objective
1280	int numberColumns = model2->numberColumns();
1281	int numberRows = model2->numberRows();
1282	int numberColumnsOrig = this->numberColumns();
1283	const double *obj = this->objective();
1284	int nNon = 0;
1285	bool allOnes = true;
1286	for (int i = 0; i < numberColumnsOrig; i++) {
1287	if (obj[i]) {
1288	nNon++;
1289	if (fabs(obj[i]) != 1.0)
1290	allOnes = false;
1291	}
1292	}
1293	if (nNon <= 1 \|\| allOnes \|\| (options.getExtraInfo(1) > 0 && options.getSpecialOption(1) == 2)) {
1294	#ifdef COIN_DEVELOP
1295	printf("Forcing primal\n");
1296	#endif
1297	method = ClpSolve::usePrimal;
1298	#ifndef CLPSOLVE_ACTIONS
1299	tryItSave = (numberRows > 200 && numberColumns > 2000 && (numberColumns > 2 * numberRows \|\| (specialOptions_ & 1024) != 0)) ? 3 : 0;
1300	#else
1301	if (numberRows > 200 && numberColumns > 2000 && numberColumns > 2 * numberRows) {
1302	tryItSave = 3;
1303	} else {
1304	// If rhs also rubbish then maybe
1305	int numberRowsOrig = this->numberRows();
1306	const double *rowLower = this->rowLower();
1307	const double *rowUpper = this->rowUpper();
1308	double last = COIN_DBL_MAX;
1309	int nDifferent = 0;
1310	for (int i = 0; i < numberRowsOrig; i++) {
1311	double value = fabs(rowLower[i]);
1312	if (value && value != COIN_DBL_MAX) {
1313	if (value != last) {
1314	nDifferent++;
1315	last = value;
1316	}
1317	}
1318	value = fabs(rowUpper[i]);
1319	if (value && value != COIN_DBL_MAX) {
1320	if (value != last) {
1321	nDifferent++;
1322	last = value;
1323	}
1324	}
1325	}
1326	if (nDifferent < 2)
1327	tryItSave = 1;
1328	}
1329	#endif
1330	}
1331	}
1332	if (method != ClpSolve::useDual && method != ClpSolve::useBarrier
1333	&& method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe
1334	&& method != ClpSolve::useBarrierNoCross) {
1335	switch (options.getSpecialOption(1)) {
1336	case 0:
1337	doIdiot = -1;
1338	doCrash = -1;
1339	doSprint = -1;
1340	break;
1341	case 1:
1342	doIdiot = 0;
1343	doCrash = 1;
1344	if (options.getExtraInfo(1) > 0)
1345	doCrash = options.getExtraInfo(1);
1346	doSprint = 0;
1347	break;
1348	case 2:
1349	doIdiot = 1;
1350	if (options.getExtraInfo(1) > 0)
1351	doIdiot = options.getExtraInfo(1);
1352	doCrash = 0;
1353	doSprint = 0;
1354	break;
1355	case 3:
1356	doIdiot = 0;
1357	doCrash = 0;
1358	doSprint = 1;
1359	break;
1360	case 4:
1361	doIdiot = 0;
1362	doCrash = 0;
1363	doSprint = 0;
1364	break;
1365	case 5:
1366	doIdiot = 0;
1367	doCrash = -1;
1368	doSprint = -1;
1369	break;
1370	case 6:
1371	doIdiot = -1;
1372	doCrash = -1;
1373	doSprint = 0;
1374	break;
1375	case 7:
1376	doIdiot = -1;
1377	doCrash = 0;
1378	doSprint = -1;
1379	break;
1380	case 8:
1381	doIdiot = -1;
1382	doCrash = 0;
1383	doSprint = 0;
1384	break;
1385	case 9:
1386	doIdiot = 0;
1387	doCrash = 0;
1388	doSprint = -1;
1389	break;
1390	case 10:
1391	doIdiot = 0;
1392	doCrash = 0;
1393	doSprint = 0;
1394	if (options.getExtraInfo(1))
1395	doSlp = options.getExtraInfo(1);
1396	break;
1397	case 11:
1398	doIdiot = 0;
1399	doCrash = 0;
1400	doSprint = 0;
1401	primalStartup = 0;
1402	break;
1403	default:
1404	abort();
1405	}
1406	} else if (method != ClpSolve::tryBenders && method != ClpSolve::tryDantzigWolfe) {
1407	// Dual
1408	switch (options.getSpecialOption(0)) {
1409	case 0:
1410	doIdiot = 0;
1411	doCrash = 0;
1412	doSprint = 0;
1413	break;
1414	case 1:
1415	doIdiot = 0;
1416	doCrash = 1;
1417	if (options.getExtraInfo(0) > 0)
1418	doCrash = options.getExtraInfo(0);
1419	doSprint = 0;
1420	break;
1421	case 2:
1422	doIdiot = -1;
1423	if (options.getExtraInfo(0) > 0)
1424	doIdiot = options.getExtraInfo(0);
1425	doCrash = 0;
1426	doSprint = 0;
1427	break;
1428	default:
1429	abort();
1430	}
1431	} else {
1432	// decomposition
1433	}
1434	#ifndef NO_RTTI
1435	ClpQuadraticObjective quadraticObj = (dynamic_cast< ClpQuadraticObjective >(objectiveAsObject()));
1436	#else
1437	ClpQuadraticObjective *quadraticObj = NULL;
1438	if (objective_->type() == 2)
1439	quadraticObj = (static_cast< ClpQuadraticObjective * >(objective_));
1440	#endif
1441	// If quadratic then primal or barrier or slp
1442	if (quadraticObj) {
1443	doSprint = 0;
1444	//doIdiot = 0;
1445	// off
1446	if (method == ClpSolve::useBarrier)
1447	method = ClpSolve::useBarrierNoCross;
1448	else if (method != ClpSolve::useBarrierNoCross)
1449	method = ClpSolve::usePrimal;
1450	}
1451	#ifdef COIN_HAS_VOL
1452	// Save number of idiot
1453	int saveDoIdiot = doIdiot;
1454	#endif
1455	// Just do this number of passes in Sprint
1456	int maxSprintPass = 100;
1457	// See if worth trying +- one matrix
1458	bool plusMinus = false;
1459	CoinBigIndex numberElements = model2->getNumElements();
1460	#ifndef SLIM_CLP
1461	#ifndef NO_RTTI
1462	if (dynamic_cast< ClpNetworkMatrix * >(matrix_)) {
1463	// network - switch off stuff
1464	doIdiot = 0;
1465	if (doSprint < 0)
1466	doSprint = 0;
1467	}
1468	#else
1469	if (matrix_->type() == 11) {
1470	// network - switch off stuff
1471	doIdiot = 0;
1472	//doSprint=0;
1473	}
1474	#endif
1475	#endif
1476	int numberColumns = model2->numberColumns();
1477	int numberRows = model2->numberRows();
1478	// If not all slack basis - switch off all except sprint
1479	int numberRowsBasic = 0;
1480	int iRow;
1481	for (iRow = 0; iRow < numberRows; iRow++)
1482	if (model2->getRowStatus(iRow) == basic)
1483	numberRowsBasic++;
1484	if (numberRowsBasic < numberRows && objective_->type() < 2) {
1485	doIdiot = 0;
1486	doCrash = 0;
1487	//doSprint=0;
1488	}
1489	if (options.getSpecialOption(3) == 0) {
1490	if (numberElements > 100000)
1491	plusMinus = true;
1492	if (numberElements > 10000 && (doIdiot \|\| doSprint))
1493	plusMinus = true;
1494	} else if ((specialOptions_ & 1024) != 0) {
1495	plusMinus = true;
1496	}
1497	#ifndef SLIM_CLP
1498	// Statistics (+1,-1, other) - used to decide on strategy if not +-1
1499	CoinBigIndex statistics[3] = { -1, 0, 0 };
1500	if (plusMinus) {
1501	saveMatrix = model2->clpMatrix();
1502	#ifndef NO_RTTI
1503	ClpPackedMatrix clpMatrix = dynamic_cast< ClpPackedMatrix >(saveMatrix);
1504	#else
1505	ClpPackedMatrix *clpMatrix = NULL;
1506	if (saveMatrix->type() == 1)
1507	clpMatrix = static_cast< ClpPackedMatrix * >(saveMatrix);
1508	#endif
1509	if (clpMatrix) {
1510	ClpPlusMinusOneMatrix newMatrix = new ClpPlusMinusOneMatrix((clpMatrix->matrix()));
1511	if (newMatrix->getIndices()) {
1512	// CHECKME This test of specialOptions and the one above
1513	// don't seem compatible.
1514	#ifndef ABC_INHERIT
1515	if ((specialOptions_ & 1024) == 0) {
1516	model2->replaceMatrix(newMatrix);
1517	} else {
1518	#endif
1519	// in integer (or abc) - just use for sprint/idiot
1520	saveMatrix = NULL;
1521	delete newMatrix;
1522	#ifndef ABC_INHERIT
1523	}
1524	#endif
1525	} else {
1526	handler_->message(CLP_MATRIX_CHANGE, messages_)
1527	<< "+- 1"
1528	<< CoinMessageEol;
1529	CoinMemcpyN(newMatrix->startPositive(), 3, statistics);
1530	saveMatrix = NULL;
1531	plusMinus = false;
1532	delete newMatrix;
1533	}
1534	} else {
1535	saveMatrix = NULL;
1536	plusMinus = false;
1537	}
1538	}
1539	#endif
1540	if (this->factorizationFrequency() == 200) {
1541	// User did not touch preset
1542	model2->defaultFactorizationFrequency();
1543	} else if (model2 != this) {
1544	// make sure model2 has correct value
1545	model2->setFactorizationFrequency(this->factorizationFrequency());
1546	}
1547	if (method == ClpSolve::automatic) {
1548	if (doSprint == 0 && doIdiot == 0) {
1549	// off
1550	method = ClpSolve::useDual;
1551	} else {
1552	// only do primal if sprint or idiot
1553	if (doSprint > 0) {
1554	method = ClpSolve::usePrimalorSprint;
1555	} else if (doIdiot > 0) {
1556	method = ClpSolve::usePrimal;
1557	} else {
1558	if (numberElements < 500000) {
1559	// Small problem
1560	if (numberRows * 10 > numberColumns \|\| numberColumns < 6000
1561	\|\| (numberRows * 20 > numberColumns && !plusMinus))
1562	doSprint = 0; // switch off sprint
1563	} else {
1564	// larger problem
1565	if (numberRows * 8 > numberColumns)
1566	doSprint = 0; // switch off sprint
1567	}
1568	// switch off idiot or sprint if any free variable
1569	// switch off sprint if very few with costs
1570	// or great variation in cost
1571	int iColumn;
1572	const double *columnLower = model2->columnLower();
1573	const double *columnUpper = model2->columnUpper();
1574	const double *objective = model2->objective();
1575	int nObj = 0;
1576	int nFree = 0;
1577	double smallestObj = COIN_DBL_MAX;
1578	double largestObj = 0.0;
1579	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1580	if (columnLower[iColumn] < -1.0e10 && columnUpper[iColumn] > 1.0e10) {
1581	nFree++;
1582	} else if (objective[iColumn]) {
1583	nObj++;
1584	smallestObj = CoinMin(smallestObj, objective[iColumn]);
1585	largestObj = CoinMax(largestObj, objective[iColumn]);
1586	}
1587	}
1588	if (nObj * 10 < numberColumns \|\| smallestObj * 10.0 < largestObj)
1589	doSprint = 0;
1590	if (nFree)
1591	doIdiot = 0;
1592	if (nFree * 10 > numberRows)
1593	doSprint = 0;
1594	int nPasses = 0;
1595	// look at rhs
1596	int iRow;
1597	double largest = 0.0;
1598	double smallest = 1.0e30;
1599	double largestGap = 0.0;
1600	int numberNotE = 0;
1601	bool notInteger = false;
1602	for (iRow = 0; iRow < numberRows; iRow++) {
1603	double value1 = model2->rowLower_[iRow];
1604	if (value1 && value1 > -1.0e31) {
1605	largest = CoinMax(largest, fabs(value1));
1606	smallest = CoinMin(smallest, fabs(value1));
1607	if (fabs(value1 - floor(value1 + 0.5)) > 1.0e-8) {
1608	notInteger = true;
1609	break;
1610	}
1611	}
1612	double value2 = model2->rowUpper_[iRow];
1613	if (value2 && value2 < 1.0e31) {
1614	largest = CoinMax(largest, fabs(value2));
1615	smallest = CoinMin(smallest, fabs(value2));
1616	if (fabs(value2 - floor(value2 + 0.5)) > 1.0e-8) {
1617	notInteger = true;
1618	break;
1619	}
1620	}
1621	// CHECKME This next bit can't be right...
1622	if (value2 > value1) {
1623	numberNotE++;
1624	//if (value2 > 1.0e31 \|\| value1 < -1.0e31)
1625	// largestGap = COIN_DBL_MAX;
1626	//else
1627	// largestGap = value2 - value1;
1628	}
1629	}
1630	int tryIt = (numberRows > 200 && numberColumns > 2000 && (numberColumns > 2 * numberRows \|\| (method != ClpSolve::useDual && (specialOptions_ & 1024) != 0))) ? 3 : 0;
1631	tryItSave = tryIt;
1632	if (numberRows < 1000 && numberColumns < 3000)
1633	tryIt = 0;
1634	if (notInteger)
1635	tryIt = 0;
1636	if (largest / smallest > 10 \|\| (largest / smallest > 2.0 && largest > 50))
1637	tryIt = 0;
1638	if (tryIt) {
1639	if (largest / smallest > 2.0) {
1640	nPasses = 10 + numberColumns / 100000;
1641	nPasses = CoinMin(nPasses, 50);
1642	nPasses = CoinMax(nPasses, 15);
1643	if (numberRows > 20000 && nPasses > 5) {
1644	// Might as well go for it
1645	nPasses = CoinMax(nPasses, 71);
1646	} else if (numberRows > 2000 && nPasses > 5) {
1647	nPasses = CoinMax(nPasses, 50);
1648	} else if (numberElements < 3 * numberColumns) {
1649	nPasses = CoinMin(nPasses, 10); // probably not worh it
1650	}
1651	} else if (largest / smallest > 1.01 \|\| numberElements <= 3 * numberColumns) {
1652	nPasses = 10 + numberColumns / 1000;
1653	nPasses = CoinMin(nPasses, 100);
1654	nPasses = CoinMax(nPasses, 30);
1655	if (numberRows > 25000) {
1656	// Might as well go for it
1657	nPasses = CoinMax(nPasses, 71);
1658	}
1659	if (!largestGap)
1660	nPasses *= 2;
1661	} else {
1662	nPasses = 10 + numberColumns / 1000;
1663	nPasses = CoinMax(nPasses, 100);
1664	if (!largestGap)
1665	nPasses *= 2;
1666	nPasses = CoinMin(nPasses, 200);
1667	}
1668	}
1669	//printf("%d rows %d cols plus %c tryIt %c largest %g smallest %g largestGap %g npasses %d sprint %c\n",
1670	// numberRows,numberColumns,plusMinus ? 'Y' : 'N',
1671	// tryIt ? 'Y' :'N',largest,smallest,largestGap,nPasses,doSprint ? 'Y' :'N');
1672	//exit(0);
1673	if (!tryIt \|\| nPasses <= 5)
1674	doIdiot = 0;
1675	#if CLPSOLVE_ACTIONS
1676	if (doIdiot && doSprint < 0 && wasAutomatic && 20 * model2->numberRows() > model2->numberColumns())
1677	doSprint = 0; // switch off sprint
1678	#endif
1679	if (doSprint) {
1680	method = ClpSolve::usePrimalorSprint;
1681	} else if (doIdiot) {
1682	method = ClpSolve::usePrimal;
1683	} else {
1684	method = ClpSolve::useDual;
1685	}
1686	}
1687	}
1688	}
1689	if (method == ClpSolve::tryBenders) {
1690	// Now build model
1691	int lengthNames = model2->lengthNames();
1692	model2->setLengthNames(0);
1693	CoinModel *build = model2->createCoinModel();
1694	model2->setLengthNames(lengthNames);
1695	CoinStructuredModel benders;
1696	build->convertMatrix();
1697	int numberBlocks = options.independentOption(0);
1698	benders.setMessageHandler(handler_);
1699	numberBlocks = benders.decompose(*build, 2, numberBlocks, NULL);
1700	delete build;
1701	//exit(0);
1702	if (numberBlocks) {
1703	options.setIndependentOption(1, 1); // don't do final clean up
1704	model2->solveBenders(&benders, options);
1705	//move solution
1706	method = ClpSolve::notImplemented;
1707	time2 = CoinCpuTime();
1708	timeCore = time2 - timeX;
1709	handler_->message(CLP_INTERVAL_TIMING, messages_)
1710	<< "Crossover" << timeCore << time2 - time1
1711	<< CoinMessageEol;
1712	timeX = time2;
1713	} else {
1714	printf("No structure\n");
1715	method = ClpSolve::useDual;
1716	}
1717	} else if (method == ClpSolve::tryDantzigWolfe) {
1718	abort();
1719	}
1720	if (method == ClpSolve::usePrimalorSprint) {
1721	if (doSprint < 0) {
1722	if (numberElements < 500000) {
1723	// Small problem
1724	if (numberRows * 10 > numberColumns \|\| numberColumns < 6000
1725	\|\| (numberRows * 20 > numberColumns && !plusMinus))
1726	method = ClpSolve::usePrimal; // switch off sprint
1727	} else {
1728	// larger problem
1729	if (numberRows * 8 > numberColumns)
1730	method = ClpSolve::usePrimal; // switch off sprint
1731	// but make lightweight
1732	if (numberRows * 10 > numberColumns \|\| numberColumns < 6000
1733	\|\| (numberRows * 20 > numberColumns && !plusMinus))
1734	maxSprintPass = 10;
1735	}
1736	} else if (doSprint == 0) {
1737	method = ClpSolve::usePrimal; // switch off sprint
1738	}
1739	}
1740	if (method == ClpSolve::useDual) {
1741	#ifdef CLP_USEFUL_PRINTOUT
1742	debugInt[6] = 1;
1743	#endif
1744	double *saveLower = NULL;
1745	double *saveUpper = NULL;
1746	if (presolve == ClpSolve::presolveOn) {
1747	int numberInfeasibilities = model2->tightenPrimalBounds(0.0, 0);
1748	if (numberInfeasibilities) {
1749	handler_->message(CLP_INFEASIBLE, messages_)
1750	<< CoinMessageEol;
1751	delete model2;
1752	model2 = this;
1753	presolve = ClpSolve::presolveOff;
1754	}
1755	} else if (numberRows_ + numberColumns_ > 5000) {
1756	// do anyway
1757	saveLower = new double[numberRows_ + numberColumns_];
1758	CoinMemcpyN(model2->columnLower(), numberColumns_, saveLower);
1759	CoinMemcpyN(model2->rowLower(), numberRows_, saveLower + numberColumns_);
1760	saveUpper = new double[numberRows_ + numberColumns_];
1761	CoinMemcpyN(model2->columnUpper(), numberColumns_, saveUpper);
1762	CoinMemcpyN(model2->rowUpper(), numberRows_, saveUpper + numberColumns_);
1763	int numberInfeasibilities = model2->tightenPrimalBounds();
1764	if (numberInfeasibilities) {
1765	handler_->message(CLP_INFEASIBLE, messages_)
1766	<< CoinMessageEol;
1767	CoinMemcpyN(saveLower, numberColumns_, model2->columnLower());
1768	CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower());
1769	delete[] saveLower;
1770	saveLower = NULL;
1771	CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper());
1772	CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper());
1773	delete[] saveUpper;
1774	saveUpper = NULL;
1775	// return if wanted
1776	if (options.infeasibleReturn() \|\| (moreSpecialOptions_ & 1) != 0)
1777	return -1;
1778	}
1779	}
1780	#ifndef COIN_HAS_VOL
1781	// switch off idiot and volume for now
1782	doIdiot = 0;
1783	#endif
1784	// pick up number passes
1785	int nPasses = 0;
1786	int numberNotE = 0;
1787	#ifndef SLIM_CLP
1788	if ((doIdiot < 0 && plusMinus) \|\| doIdiot > 0) {
1789	// See if candidate for idiot
1790	nPasses = 0;
1791	Idiot info(*model2);
1792	info.setStrategy(idiotOptions \| info.getStrategy());
1793	// Get average number of elements per column
1794	double ratio = static_cast< double >(numberElements) / static_cast< double >(numberColumns);
1795	// look at rhs
1796	int iRow;
1797	double largest = 0.0;
1798	double smallest = 1.0e30;
1799	for (iRow = 0; iRow < numberRows; iRow++) {
1800	double value1 = model2->rowLower_[iRow];
1801	if (value1 && value1 > -1.0e31) {
1802	largest = CoinMax(largest, fabs(value1));
1803	smallest = CoinMin(smallest, fabs(value1));
1804	}
1805	double value2 = model2->rowUpper_[iRow];
1806	if (value2 && value2 < 1.0e31) {
1807	largest = CoinMax(largest, fabs(value2));
1808	smallest = CoinMin(smallest, fabs(value2));
1809	}
1810	if (value2 > value1) {
1811	numberNotE++;
1812	}
1813	}
1814	if (doIdiot < 0) {
1815	if (numberRows > 200 && numberColumns > 5000 && ratio >= 3.0 && largest / smallest < 1.1 && !numberNotE) {
1816	nPasses = 71;
1817	}
1818	}
1819	if (doIdiot > 0) {
1820	nPasses = CoinMax(nPasses, doIdiot);
1821	if (nPasses > 70) {
1822	info.setStartingWeight(1.0e3);
1823	info.setDropEnoughFeasibility(0.01);
1824	}
1825	}
1826	if (nPasses > 20) {
1827	#ifdef COIN_HAS_VOL
1828	int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot);
1829	if (!returnCode) {
1830	time2 = CoinCpuTime();
1831	timeIdiot = time2 - timeX;
1832	handler_->message(CLP_INTERVAL_TIMING, messages_)
1833	<< "Idiot Crash" << timeIdiot << time2 - time1
1834	<< CoinMessageEol;
1835	timeX = time2;
1836	} else {
1837	nPasses = 0;
1838	}
1839	#else
1840	nPasses = 0;
1841	#endif
1842	} else {
1843	nPasses = 0;
1844	}
1845	}
1846	#endif
1847	if (doCrash) {
1848	#ifdef ABC_INHERIT
1849	if (!model2->abcState()) {
1850	#endif
1851	switch (doCrash) {
1852	// standard
1853	case 1:
1854	model2->crash(1000, 1);
1855	break;
1856	// As in paper by Solow and Halim (approx)
1857	case 2:
1858	case 3:
1859	model2->crash(model2->dualBound(), 0);
1860	break;
1861	// Just put free in basis
1862	case 4:
1863	model2->crash(0.0, 3);
1864	break;
1865	}
1866	#ifdef ABC_INHERIT
1867	} else if (doCrash >= 0) {
1868	model2->setAbcState(model2->abcState() \| 256 * doCrash);
1869	}
1870	#endif
1871	}
1872	if (!nPasses) {
1873	int saveOptions = model2->specialOptions();
1874	if (model2->numberRows() > 100)
1875	model2->setSpecialOptions(saveOptions \| 64); // go as far as possible
1876	//int numberRows = model2->numberRows();
1877	//int numberColumns = model2->numberColumns();
1878	if (dynamic_cast< ClpPackedMatrix * >(matrix_)) {
1879	// See if original wanted vector
1880	ClpPackedMatrix clpMatrixO = dynamic_cast< ClpPackedMatrix >(matrix_);
1881	ClpMatrixBase *matrix = model2->clpMatrix();
1882	if (dynamic_cast< ClpPackedMatrix * >(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
1883	ClpPackedMatrix clpMatrix = dynamic_cast< ClpPackedMatrix >(matrix);
1884	clpMatrix->makeSpecialColumnCopy();
1885	//model2->setSpecialOptions(model2->specialOptions()\|256); // to say no row copy for comparisons
1886	model2->dual(0);
1887	clpMatrix->releaseSpecialColumnCopy();
1888	} else {
1889	#ifndef ABC_INHERIT
1890	model2->dual(0);
1891	#else
1892	model2->dealWithAbc(0, 0, interrupt);
1893	#endif
1894	}
1895	} else {
1896	model2->dual(0);
1897	}
1898	} else if (!numberNotE && 0) {
1899	// E so we can do in another way
1900	double *pi = model2->dualRowSolution();
1901	int i;
1902	int numberColumns = model2->numberColumns();
1903	int numberRows = model2->numberRows();
1904	double *saveObj = new double[numberColumns];
1905	CoinMemcpyN(model2->objective(), numberColumns, saveObj);
1906	CoinMemcpyN(model2->objective(),
1907	numberColumns, model2->dualColumnSolution());
1908	model2->clpMatrix()->transposeTimes(-1.0, pi, model2->dualColumnSolution());
1909	CoinMemcpyN(model2->dualColumnSolution(),
1910	numberColumns, model2->objective());
1911	const double *rowsol = model2->primalRowSolution();
1912	double offset = 0.0;
1913	for (i = 0; i < numberRows; i++) {
1914	offset += pi[i] * rowsol[i];
1915	}
1916	double value2;
1917	model2->getDblParam(ClpObjOffset, value2);
1918	//printf("Offset %g %g\n",offset,value2);
1919	model2->setDblParam(ClpObjOffset, value2 - offset);
1920	model2->setPerturbation(51);
1921	//model2->setRowObjective(pi);
1922	// zero out pi
1923	//memset(pi,0,numberRows*sizeof(double));
1924	// Could put some in basis - only partially tested
1925	model2->allSlackBasis();
1926	//model2->factorization()->maximumPivots(200);
1927	//model2->setLogLevel(63);
1928	// solve
1929	model2->dual(0);
1930	model2->setDblParam(ClpObjOffset, value2);
1931	CoinMemcpyN(saveObj, numberColumns, model2->objective());
1932	// zero out pi
1933	//memset(pi,0,numberRows*sizeof(double));
1934	//model2->setRowObjective(pi);
1935	delete[] saveObj;
1936	//model2->dual(0);
1937	model2->setPerturbation(50);
1938	model2->primal();
1939	} else {
1940	// solve
1941	model2->setPerturbation(100);
1942	model2->dual(2);
1943	model2->setPerturbation(50);
1944	model2->dual(0);
1945	}
1946	if (saveLower) {
1947	CoinMemcpyN(saveLower, numberColumns_, model2->columnLower());
1948	CoinMemcpyN(saveLower + numberColumns_, numberRows_, model2->rowLower());
1949	delete[] saveLower;
1950	saveLower = NULL;
1951	CoinMemcpyN(saveUpper, numberColumns_, model2->columnUpper());
1952	CoinMemcpyN(saveUpper + numberColumns_, numberRows_, model2->rowUpper());
1953	delete[] saveUpper;
1954	saveUpper = NULL;
1955	}
1956	time2 = CoinCpuTime();
1957	timeCore = time2 - timeX;
1958	handler_->message(CLP_INTERVAL_TIMING, messages_)
1959	<< "Dual" << timeCore << time2 - time1
1960	<< CoinMessageEol;
1961	timeX = time2;
1962	} else if (method == ClpSolve::usePrimal) {
1963	#ifdef CLP_USEFUL_PRINTOUT
1964	debugInt[6] = 2;
1965	#endif
1966	#ifndef SLIM_CLP
1967	if (doIdiot) {
1968	int nPasses = 0;
1969	Idiot info(*model2);
1970	info.setStrategy(idiotOptions \| info.getStrategy());
1971	// Get average number of elements per column
1972	double ratio = static_cast< double >(numberElements) / static_cast< double >(numberColumns);
1973	// look at rhs
1974	int iRow;
1975	double largest = 0.0;
1976	double smallest = 1.0e30;
1977	double largestGap = 0.0;
1978	int numberNotE = 0;
1979	for (iRow = 0; iRow < numberRows; iRow++) {
1980	double value1 = model2->rowLower_[iRow];
1981	if (value1 && value1 > -1.0e31) {
1982	largest = CoinMax(largest, fabs(value1));
1983	smallest = CoinMin(smallest, fabs(value1));
1984	}
1985	double value2 = model2->rowUpper_[iRow];
1986	if (value2 && value2 < 1.0e31) {
1987	largest = CoinMax(largest, fabs(value2));
1988	smallest = CoinMin(smallest, fabs(value2));
1989	}
1990	if (value2 > value1) {
1991	numberNotE++;
1992	if (value2 > 1.0e31 \|\| value1 < -1.0e31)
1993	largestGap = COIN_DBL_MAX;
1994	else
1995	largestGap = value2 - value1;
1996	}
1997	}
1998	bool increaseSprint = plusMinus;
1999	if ((specialOptions_ & 1024) != 0)
2000	increaseSprint = false;
2001	if (!plusMinus) {
2002	// If 90% +- 1 then go for sprint
2003	if (statistics[0] >= 0 && 10 * statistics[2] < statistics[0] + statistics[1])
2004	increaseSprint = true;
2005	}
2006	int tryIt = tryItSave;
2007	if (numberRows < 1000 && numberColumns < 3000)
2008	tryIt = 0;
2009	if (tryIt) {
2010	if (increaseSprint) {
2011	info.setStartingWeight(1.0e3);
2012	info.setReduceIterations(6);
2013	// also be more lenient on infeasibilities
2014	info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility());
2015	info.setDropEnoughWeighted(-2.0);
2016	if (largest / smallest > 2.0) {
2017	nPasses = 10 + numberColumns / 100000;
2018	nPasses = CoinMin(nPasses, 50);
2019	nPasses = CoinMax(nPasses, 15);
2020	if (numberRows > 20000 && nPasses > 5) {
2021	// Might as well go for it
2022	nPasses = CoinMax(nPasses, 71);
2023	} else if (numberRows > 2000 && nPasses > 5) {
2024	nPasses = CoinMax(nPasses, 50);
2025	} else if (numberElements < 3 * numberColumns) {
2026	nPasses = CoinMin(nPasses, 10); // probably not worh it
2027	if (doIdiot < 0)
2028	info.setLightweight(1); // say lightweight idiot
2029	} else {
2030	if (doIdiot < 0)
2031	info.setLightweight(1); // say lightweight idiot
2032	}
2033	} else if (largest / smallest > 1.01 \|\| numberElements <= 3 * numberColumns) {
2034	nPasses = 10 + numberColumns / 1000;
2035	nPasses = CoinMin(nPasses, 100);
2036	nPasses = CoinMax(nPasses, 30);
2037	if (numberRows > 25000) {
2038	// Might as well go for it
2039	nPasses = CoinMax(nPasses, 71);
2040	}
2041	if (!largestGap)
2042	nPasses *= 2;
2043	} else {
2044	nPasses = 10 + numberColumns / 1000;
2045	nPasses = CoinMin(nPasses, 200);
2046	nPasses = CoinMax(nPasses, 100);
2047	info.setStartingWeight(1.0e-1);
2048	info.setReduceIterations(6);
2049	if (!largestGap && nPasses <= 50)
2050	nPasses *= 2;
2051	//info.setFeasibilityTolerance(1.0e-7);
2052	}
2053	// If few passes - don't bother
2054	if (nPasses <= 5 && !plusMinus)
2055	nPasses = 0;
2056	} else {
2057	if (doIdiot < 0)
2058	info.setLightweight(1); // say lightweight idiot
2059	if (largest / smallest > 1.01 \|\| numberNotE \|\| statistics[2] > statistics[0] + statistics[1]) {
2060	if (numberRows > 25000 \|\| numberColumns > 5 * numberRows) {
2061	nPasses = 50;
2062	} else if (numberColumns > 4 * numberRows) {
2063	nPasses = 20;
2064	} else {
2065	nPasses = 5;
2066	}
2067	} else {
2068	if (numberRows > 25000 \|\| numberColumns > 5 * numberRows) {
2069	nPasses = 50;
2070	info.setLightweight(0); // say not lightweight idiot
2071	} else if (numberColumns > 4 * numberRows) {
2072	nPasses = 20;
2073	} else {
2074	nPasses = 15;
2075	}
2076	}
2077	if (ratio < 3.0) {
2078	nPasses = static_cast< int >(ratio * static_cast< double >(nPasses) / 4.0); // probably not worth it
2079	} else {
2080	nPasses = CoinMax(nPasses, 5);
2081	}
2082	if (numberRows > 25000 && nPasses > 5) {
2083	// Might as well go for it
2084	nPasses = CoinMax(nPasses, 71);
2085	} else if (increaseSprint) {
2086	nPasses *= 2;
2087	nPasses = CoinMin(nPasses, 71);
2088	} else if (nPasses == 5 && ratio > 5.0) {
2089	nPasses = static_cast< int >(static_cast< double >(nPasses) * (ratio / 5.0)); // increase if lots of elements per column
2090	}
2091	if (nPasses <= 5 && !plusMinus)
2092	nPasses = 0;
2093	//info.setStartingWeight(1.0e-1);
2094	}
2095	if (tryIt == 1) {
2096	idiotOptions \|= 262144;
2097	info.setStrategy(idiotOptions \| info.getStrategy());
2098	//model2->setSpecialOptions(model2->specialOptions()
2099	// \|8388608);
2100	}
2101	}
2102	if (doIdiot > 0) {
2103	// pick up number passes
2104	nPasses = options.getExtraInfo(1) % 1000000;
2105	#ifdef COIN_HAS_VOL
2106	int returnCode = solveWithVolume(model2, nPasses, saveDoIdiot);
2107	nPasses = 0;
2108	if (!returnCode) {
2109	time2 = CoinCpuTime();
2110	timeIdiot = time2 - timeX;
2111	handler_->message(CLP_INTERVAL_TIMING, messages_)
2112	<< "Idiot Crash" << timeIdiot << time2 - time1
2113	<< CoinMessageEol;
2114	timeX = time2;
2115	}
2116	#endif
2117	#ifdef CLP_USEFUL_PRINTOUT
2118	debugInt[6] = 4;
2119	debugInt[7] = nPasses;
2120	#endif
2121	if (nPasses > 70) {
2122	info.setStartingWeight(1.0e3);
2123	info.setReduceIterations(6);
2124	//if (nPasses > 200)
2125	//info.setFeasibilityTolerance(1.0e-9);
2126	//if (nPasses > 1900)
2127	//info.setWeightFactor(0.93);
2128	if (nPasses > 900) {
2129	double reductions = nPasses / 6.0;
2130	if (nPasses < 5000) {
2131	reductions /= 12.0;
2132	} else {
2133	reductions /= 13.0;
2134	info.setStartingWeight(1.0e4);
2135	}
2136	double ratio = 1.0 / std::pow(10.0, (1.0 / reductions));
2137	printf("%d passes reduction factor %g\n", nPasses, ratio);
2138	info.setWeightFactor(ratio);
2139	} else if (nPasses > 500) {
2140	info.setWeightFactor(0.7);
2141	} else if (nPasses > 200) {
2142	info.setWeightFactor(0.5);
2143	}
2144	if (maximumIterations() < nPasses) {
2145	printf("Presuming maximumIterations is just for Idiot\n");
2146	nPasses = maximumIterations();
2147	setMaximumIterations(COIN_INT_MAX);
2148	model2->setMaximumIterations(COIN_INT_MAX);
2149	}
2150	if (nPasses >= 10000 && nPasses < 100000) {
2151	int k = nPasses % 100;
2152	nPasses /= 200;
2153	info.setReduceIterations(3);
2154	if (k)
2155	info.setStartingWeight(1.0e2);
2156	}
2157	// also be more lenient on infeasibilities
2158	info.setDropEnoughFeasibility(0.5 * info.getDropEnoughFeasibility());
2159	info.setDropEnoughWeighted(-2.0);
2160	} else if (nPasses >= 50) {
2161	info.setStartingWeight(1.0e3);
2162	//info.setReduceIterations(6);
2163	}
2164	// For experimenting
2165	if (nPasses < 70 && (nPasses % 10) > 0 && (nPasses % 10) < 4) {
2166	info.setStartingWeight(1.0e3);
2167	info.setLightweight(nPasses % 10); // special testing
2168	#ifdef COIN_DEVELOP
2169	printf("warning - odd lightweight %d\n", nPasses % 10);
2170	//info.setReduceIterations(6);
2171	#endif
2172	}
2173	}
2174	if (options.getExtraInfo(1) > 1000000)
2175	nPasses += 1000000;
2176	if (nPasses) {
2177	doCrash = 0;
2178	#if 0
2179	double * solution = model2->primalColumnSolution();
2180	int iColumn;
2181	double * saveLower = new double[numberColumns];
2182	CoinMemcpyN(model2->columnLower(), numberColumns, saveLower);
2183	double * saveUpper = new double[numberColumns];
2184	CoinMemcpyN(model2->columnUpper(), numberColumns, saveUpper);
2185	printf("doing tighten before idiot\n");
2186	model2->tightenPrimalBounds();
2187	// Move solution
2188	double * columnLower = model2->columnLower();
2189	double * columnUpper = model2->columnUpper();
2190	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2191	if (columnLower[iColumn] > 0.0)
2192	solution[iColumn] = columnLower[iColumn];
2193	else if (columnUpper[iColumn] < 0.0)
2194	solution[iColumn] = columnUpper[iColumn];
2195	else
2196	solution[iColumn] = 0.0;
2197	}
2198	CoinMemcpyN(saveLower, numberColumns, columnLower);
2199	CoinMemcpyN(saveUpper, numberColumns, columnUpper);
2200	delete [] saveLower;
2201	delete [] saveUpper;
2202	#else
2203	// Allow for crossover
2204	//#define LACI_TRY
2205	#ifndef LACI_TRY
2206	//if (doIdiot>0)
2207	#if 0 //def ABC_INHERIT
2208	if (!model2->abcState())
2209	#endif
2210	info.setStrategy(512 \| info.getStrategy());
2211	#endif
2212	// Allow for scaling
2213	info.setStrategy(32 \| info.getStrategy());
2214	int saveScalingFlag = model2->scalingFlag();
2215	bool linearObjective = objective_->type() < 2;
2216	if (!linearObjective)
2217	model2->scaling(0);
2218	#define CLP_DUAL_IDIOT
2219	#ifdef CLP_DUAL_IDIOT
2220	bool doubleIdiot = false;
2221	if (nPasses == 99 && linearObjective)
2222	doubleIdiot = true;
2223	if (doubleIdiot) {
2224	ClpSimplex dualModel2 = static_cast< ClpSimplexOther >(model2)->dualOfModel(1.0, 1.0);
2225	if (dualModel2) {
2226	printf("Dual of model has %d rows and %d columns\n",
2227	dualModel2->numberRows(), dualModel2->numberColumns());
2228	dualModel2->setOptimizationDirection(1.0);
2229	Idiot infoDual(info);
2230	infoDual.setModel(dualModel2);
2231	#if 0
2232	info.setStrategy(512 \| info.getStrategy());
2233	// Allow for scaling
2234	info.setStrategy(32 \| info.getStrategy());
2235	info.setStartingWeight(1.0e3);
2236	info.setReduceIterations(6);
2237	#endif
2238	info.crash(nPasses, model2->messageHandler(),
2239	model2->messagesPointer(), false);
2240	infoDual.crash(nPasses, model2->messageHandler(),
2241	model2->messagesPointer(), false);
2242	// two copies of solutions
2243	ClpSimplex temp(*model2);
2244	#if 0
2245	static_cast<ClpSimplexOther *> (&temp)->restoreFromDual(dualModel2);
2246	#else
2247	// move duals and just copy primal
2248	memcpy(temp.dualRowSolution(), dualModel2->primalColumnSolution(),
2249	numberRows * sizeof(double));
2250	memcpy(temp.primalColumnSolution(), model2->primalColumnSolution(),
2251	numberColumns * sizeof(double));
2252	#endif
2253	delete dualModel2;
2254	int numberRows = model2->numberRows();
2255	int numberColumns = model2->numberColumns();
2256	ClpSimplex *tempModel[2];
2257	tempModel[0] = model2;
2258	tempModel[1] = &temp;
2259	const double *primalColumn[2];
2260	const double *dualRow[2];
2261	double *dualColumn[2];
2262	double *primalRow[2];
2263	for (int i = 0; i < 2; i++) {
2264	primalColumn[i] = tempModel[i]->primalColumnSolution();
2265	dualRow[i] = tempModel[i]->dualRowSolution();
2266	dualColumn[i] = tempModel[i]->dualColumnSolution();
2267	primalRow[i] = tempModel[i]->primalRowSolution();
2268	memcpy(dualColumn[i], model2->objective(),
2269	numberColumns * sizeof(double));
2270	memset(primalRow[i], 0, numberRows * sizeof(double));
2271	tempModel[i]->clpMatrix()->transposeTimes(-1.0,
2272	dualRow[i],
2273	dualColumn[i]);
2274	tempModel[i]->clpMatrix()->times(1.0,
2275	primalColumn[i],
2276	primalRow[i]);
2277	tempModel[i]->checkSolutionInternal();
2278	printf("model %d - dual inf %g primal inf %g\n",
2279	i, tempModel[i]->sumDualInfeasibilities(),
2280	tempModel[i]->sumPrimalInfeasibilities());
2281	}
2282	printf("What now\n");
2283	} else {
2284	doubleIdiot = false;
2285	}
2286	}
2287	if (!doubleIdiot)
2288	info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(),
2289	(objective_->type() < 2));
2290	#else
2291	info.crash(nPasses, model2->messageHandler(), model2->messagesPointer(), (objective_->type() < 2));
2292	#endif
2293	model2->scaling(saveScalingFlag);
2294	#endif
2295	time2 = CoinCpuTime();
2296	timeIdiot = time2 - timeX;
2297	handler_->message(CLP_INTERVAL_TIMING, messages_)
2298	<< "Idiot Crash" << timeIdiot << time2 - time1
2299	<< CoinMessageEol;
2300	timeX = time2;
2301	if (nPasses > 100000 && nPasses < 100500) {
2302	// make sure no status left
2303	model2->createStatus();
2304	// solve
2305	if (model2->factorizationFrequency() == 200) {
2306	// User did not touch preset
2307	model2->defaultFactorizationFrequency();
2308	}
2309	//int numberRows = model2->numberRows();
2310	int numberColumns = model2->numberColumns();
2311	// save duals
2312	//double * saveDuals = CoinCopyOfArray(model2->dualRowSolution(),numberRows);
2313	// for moment this only works on nug etc (i.e. all ==)
2314	// needs row objective
2315	double *saveObj = CoinCopyOfArray(model2->objective(), numberColumns);
2316	double *pi = model2->dualRowSolution();
2317	model2->clpMatrix()->transposeTimes(-1.0, pi, model2->objective());
2318	// just primal values pass
2319	double saveScale = model2->objectiveScale();
2320	model2->setObjectiveScale(1.0e-3);
2321	model2->primal(2);
2322	model2->writeMps("xx.mps");
2323	double *solution = model2->primalColumnSolution();
2324	double *upper = model2->columnUpper();
2325	for (int i = 0; i < numberColumns; i++) {
2326	if (solution[i] < 100.0)
2327	upper[i] = 1000.0;
2328	}
2329	model2->setProblemStatus(-1);
2330	model2->setObjectiveScale(saveScale);
2331	#ifdef ABC_INHERIT
2332	AbcSimplex abcModel2 = new AbcSimplex(model2);
2333	if (interrupt)
2334	currentAbcModel = abcModel2;
2335	if (abcSimplex_) {
2336	// move factorization stuff
2337	abcModel2->factorization()->synchronize(model2->factorization(), abcModel2);
2338	}
2339	abcModel2->startPermanentArrays();
2340	abcModel2->setAbcState(CLP_ABC_WANTED);
2341	#if ABC_PARALLEL
2342	int parallelMode = 1;
2343	printf("Parallel mode %d\n", parallelMode);
2344	abcModel2->setParallelMode(parallelMode);
2345	#endif
2346	//if (processTune>0&&processTune<8)
2347	//abcModel2->setMoreSpecialOptions(abcModel2->moreSpecialOptions()\|65536*processTune);
2348	abcModel2->doAbcDual();
2349	abcModel2->moveStatusToClp(model2);
2350	//ClpModel::stopPermanentArrays();
2351	model2->setSpecialOptions(model2->specialOptions() & ~65536);
2352	//model2->dual();
2353	//model2->setNumberIterations(abcModel2->numberIterations()+model2->numberIterations());
2354	delete abcModel2;
2355	#endif
2356	memcpy(model2->objective(), saveObj, numberColumns * sizeof(double));
2357	//delete [] saveDuals;
2358	delete[] saveObj;
2359	model2->dual(2);
2360	} // end dubious idiot
2361	}
2362	}
2363	#endif
2364	// ?
2365	if (doCrash) {
2366	switch (doCrash) {
2367	// standard
2368	case 1:
2369	model2->crash(1000, 1);
2370	break;
2371	// As in paper by Solow and Halim (approx)
2372	case 2:
2373	model2->crash(model2->dualBound(), 0);
2374	break;
2375	// My take on it
2376	case 3:
2377	model2->crash(model2->dualBound(), -1);
2378	break;
2379	// Just put free in basis
2380	case 4:
2381	model2->crash(0.0, 3);
2382	break;
2383	}
2384	}
2385	#ifndef SLIM_CLP
2386	if (doSlp && objective_->type() == 2) {
2387	model2->nonlinearSLP(doSlp, 1.0e-5);
2388	}
2389	#endif
2390	#ifndef LACI_TRY
2391	if (options.getSpecialOption(1) != 2 \|\| options.getExtraInfo(1) < 1000000) {
2392	if (dynamic_cast< ClpPackedMatrix * >(matrix_)) {
2393	// See if original wanted vector
2394	ClpPackedMatrix clpMatrixO = dynamic_cast< ClpPackedMatrix >(matrix_);
2395	ClpMatrixBase *matrix = model2->clpMatrix();
2396	if (dynamic_cast< ClpPackedMatrix * >(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
2397	ClpPackedMatrix clpMatrix = dynamic_cast< ClpPackedMatrix >(matrix);
2398	clpMatrix->makeSpecialColumnCopy();
2399	//model2->setSpecialOptions(model2->specialOptions()\|256); // to say no row copy for comparisons
2400	model2->primal(primalStartup);
2401	clpMatrix->releaseSpecialColumnCopy();
2402	} else {
2403	#ifndef ABC_INHERIT
2404	model2->primal(primalStartup);
2405	#else
2406	model2->dealWithAbc(1, primalStartup, interrupt);
2407	#endif
2408	}
2409	} else {
2410	#ifndef ABC_INHERIT
2411	model2->primal(primalStartup);
2412	#else
2413	model2->dealWithAbc(1, primalStartup, interrupt);
2414	#endif
2415	}
2416	}
2417	#endif
2418	time2 = CoinCpuTime();
2419	timeCore = time2 - timeX;
2420	handler_->message(CLP_INTERVAL_TIMING, messages_)
2421	<< "Primal" << timeCore << time2 - time1
2422	<< CoinMessageEol;
2423	timeX = time2;
2424	} else if (method == ClpSolve::usePrimalorSprint) {
2425	// Sprint
2426	/*
2427	This driver implements what I called Sprint when I introduced the idea
2428	many years ago. Cplex calls it "sifting" which I think is just as silly.
2429	When I thought of this trivial idea
2430	it reminded me of an LP code of the 60's called sprint which after
2431	every factorization took a subset of the matrix into memory (all
2432	64K words!) and then iterated very fast on that subset. On the
2433	problems of those days it did not work very well, but it worked very
2434	well on aircrew scheduling problems where there were very large numbers
2435	of columns all with the same flavor.
2436	*/
2437
2438	/* The idea works best if you can get feasible easily. To make it
2439	more general we can add in costed slacks */
2440
2441	int originalNumberColumns = model2->numberColumns();
2442	int numberRows = model2->numberRows();
2443	ClpSimplex *originalModel2 = model2;
2444
2445	// We will need arrays to choose variables. These are too big but ..
2446	double *weight = new double[numberRows + originalNumberColumns];
2447	int *sort = new int[numberRows + originalNumberColumns];
2448	int numberSort = 0;
2449	// We are going to add slacks to get feasible.
2450	// initial list will just be artificials
2451	int iColumn;
2452	const double *columnLower = model2->columnLower();
2453	const double *columnUpper = model2->columnUpper();
2454	double *columnSolution = model2->primalColumnSolution();
2455
2456	// See if we have costed slacks
2457	int *negSlack = new int[numberRows + numberColumns];
2458	int *posSlack = new int[numberRows + numberColumns];
2459	int *nextNegSlack = negSlack + numberRows;
2460	int *nextPosSlack = posSlack + numberRows;
2461	int iRow;
2462	for (iRow = 0; iRow < numberRows + numberColumns; iRow++) {
2463	negSlack[iRow] = -1;
2464	posSlack[iRow] = -1;
2465	}
2466	const double *element = model2->matrix()->getElements();
2467	const int *row = model2->matrix()->getIndices();
2468	const CoinBigIndex *columnStart = model2->matrix()->getVectorStarts();
2469	const int *columnLength = model2->matrix()->getVectorLengths();
2470	//bool allSlack = (numberRowsBasic==numberRows);
2471	for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
2472	if (!columnSolution[iColumn] \|\| fabs(columnSolution[iColumn]) > 1.0e20) {
2473	double value = 0.0;
2474	if (columnLower[iColumn] > 0.0)
2475	value = columnLower[iColumn];
2476	else if (columnUpper[iColumn] < 0.0)
2477	value = columnUpper[iColumn];
2478	columnSolution[iColumn] = value;
2479	}
2480	if (columnLength[iColumn] == 1) {
2481	int jRow = row[columnStart[iColumn]];
2482	if (!columnLower[iColumn]) {
2483	if (element[columnStart[iColumn]] > 0.0) {
2484	if (posSlack[jRow] < 0) {
2485	posSlack[jRow] = iColumn;
2486	} else {
2487	int jColumn = posSlack[jRow];
2488	while (nextPosSlack[jColumn] >= 0)
2489	jColumn = nextPosSlack[jColumn];
2490	nextPosSlack[jColumn] = iColumn;
2491	}
2492	} else if (element[columnStart[iColumn]] < 0.0) {
2493	if (negSlack[jRow] < 0) {
2494	negSlack[jRow] = iColumn;
2495	} else {
2496	int jColumn = negSlack[jRow];
2497	while (nextNegSlack[jColumn] >= 0)
2498	jColumn = nextNegSlack[jColumn];
2499	nextNegSlack[jColumn] = iColumn;
2500	}
2501	}
2502	} else if (!columnUpper[iColumn] && false) { // out for testing
2503	if (element[columnStart[iColumn]] < 0.0 && posSlack[jRow] < 0)
2504	posSlack[jRow] = iColumn;
2505	else if (element[columnStart[iColumn]] > 0.0 && negSlack[jRow] < 0)
2506	negSlack[jRow] = iColumn;
2507	}
2508	}
2509	}
2510	// now see what that does to row solution
2511	const double *objective = model2->objective();
2512	double *rowSolution = model2->primalRowSolution();
2513	CoinZeroN(rowSolution, numberRows);
2514	model2->clpMatrix()->times(1.0, columnSolution, rowSolution);
2515	// See if we can adjust using costed slacks
2516	double penalty = CoinMax(1.0e5, CoinMin(infeasibilityCost_ * 0.01, 1.0e10)) * optimizationDirection_;
2517	const double *lower = model2->rowLower();
2518	const double *upper = model2->rowUpper();
2519	for (iRow = 0; iRow < numberRows; iRow++) {
2520	if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
2521	int jColumn = posSlack[iRow];
2522	if (jColumn >= 0) {
2523	// sort if more than one
2524	int nPos = 1;
2525	sort[0] = jColumn;
2526	weight[0] = objective[jColumn];
2527	while (nextPosSlack[jColumn] >= 0) {
2528	jColumn = nextPosSlack[jColumn];
2529	sort[nPos] = jColumn;
2530	weight[nPos++] = objective[jColumn];
2531	}
2532	if (nPos > 1) {
2533	CoinSort_2(weight, weight + nPos, sort);
2534	for (int i = 0; i < nPos; i++) {
2535	double difference = lower[iRow] - rowSolution[iRow];
2536	jColumn = sort[i];
2537	double elementValue = element[columnStart[jColumn]];
2538	assert(elementValue > 0.0);
2539	double value = columnSolution[jColumn];
2540	double movement = CoinMin(difference / elementValue, columnUpper[jColumn] - value);
2541	columnSolution[jColumn] += movement;
2542	rowSolution[iRow] += movement * elementValue;
2543	}
2544	continue;
2545	}
2546	if (jColumn < 0 \|\| columnSolution[jColumn])
2547	continue;
2548	double difference = lower[iRow] - rowSolution[iRow];
2549	double elementValue = element[columnStart[jColumn]];
2550	if (elementValue > 0.0) {
2551	double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
2552	columnSolution[jColumn] = movement;
2553	rowSolution[iRow] += movement * elementValue;
2554	} else {
2555	double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
2556	columnSolution[jColumn] = movement;
2557	rowSolution[iRow] += movement * elementValue;
2558	}
2559	}
2560	} else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
2561	int jColumn = negSlack[iRow];
2562	if (jColumn >= 0) {
2563	// sort if more than one
2564	int nNeg = 1;
2565	sort[0] = jColumn;
2566	weight[0] = objective[jColumn];
2567	while (nextNegSlack[jColumn] >= 0) {
2568	jColumn = nextNegSlack[jColumn];
2569	sort[nNeg] = jColumn;
2570	weight[nNeg++] = objective[jColumn];
2571	}
2572	if (nNeg > 1) {
2573	CoinSort_2(weight, weight + nNeg, sort);
2574	for (int i = 0; i < nNeg; i++) {
2575	double difference = rowSolution[iRow] - upper[iRow];
2576	jColumn = sort[i];
2577	double elementValue = element[columnStart[jColumn]];
2578	assert(elementValue < 0.0);
2579	double value = columnSolution[jColumn];
2580	double movement = CoinMin(difference / -elementValue, columnUpper[jColumn] - value);
2581	columnSolution[jColumn] += movement;
2582	rowSolution[iRow] += movement * elementValue;
2583	}
2584	continue;
2585	}
2586	if (jColumn < 0 \|\| columnSolution[jColumn])
2587	continue;
2588	double difference = upper[iRow] - rowSolution[iRow];
2589	double elementValue = element[columnStart[jColumn]];
2590	if (elementValue < 0.0) {
2591	double movement = CoinMin(difference / elementValue, columnUpper[jColumn]);
2592	columnSolution[jColumn] = movement;
2593	rowSolution[iRow] += movement * elementValue;
2594	} else {
2595	double movement = CoinMax(difference / elementValue, columnLower[jColumn]);
2596	columnSolution[jColumn] = movement;
2597	rowSolution[iRow] += movement * elementValue;
2598	}
2599	}
2600	}
2601	}
2602	delete[] negSlack;
2603	delete[] posSlack;
2604	int nRow = numberRows;
2605	bool network = false;
2606	if (dynamic_cast< ClpNetworkMatrix * >(matrix_)) {
2607	network = true;
2608	nRow *= 2;
2609	}
2610	CoinBigIndex *addStarts = new CoinBigIndex[nRow + 1];
2611	int *addRow = new int[nRow];
2612	double *addElement = new double[nRow];
2613	addStarts[0] = 0;
2614	int numberArtificials = 0;
2615	int numberAdd = 0;
2616	double *addCost = new double[numberRows];
2617	for (iRow = 0; iRow < numberRows; iRow++) {
2618	if (lower[iRow] > rowSolution[iRow] + 1.0e-8) {
2619	addRow[numberAdd] = iRow;
2620	addElement[numberAdd++] = 1.0;
2621	if (network) {
2622	addRow[numberAdd] = numberRows;
2623	addElement[numberAdd++] = -1.0;
2624	}
2625	addCost[numberArtificials] = penalty;
2626	numberArtificials++;
2627	addStarts[numberArtificials] = numberAdd;
2628	} else if (upper[iRow] < rowSolution[iRow] - 1.0e-8) {
2629	addRow[numberAdd] = iRow;
2630	addElement[numberAdd++] = -1.0;
2631	if (network) {
2632	addRow[numberAdd] = numberRows;
2633	addElement[numberAdd++] = 1.0;
2634	}
2635	addCost[numberArtificials] = penalty;
2636	numberArtificials++;
2637	addStarts[numberArtificials] = numberAdd;
2638	}
2639	}
2640	if (numberArtificials) {
2641	// need copy so as not to disturb original
2642	model2 = new ClpSimplex(*model2);
2643	if (network) {
2644	// network - add a null row
2645	model2->addRow(0, NULL, NULL, -COIN_DBL_MAX, COIN_DBL_MAX);
2646	numberRows++;
2647	}
2648	model2->addColumns(numberArtificials, NULL, NULL, addCost,
2649	addStarts, addRow, addElement);
2650	}
2651	delete[] addStarts;
2652	delete[] addRow;
2653	delete[] addElement;
2654	delete[] addCost;
2655	// look at rhs to see if to perturb
2656	double largest = 0.0;
2657	double smallest = 1.0e30;
2658	for (iRow = 0; iRow < numberRows; iRow++) {
2659	double value;
2660	value = fabs(model2->rowLower_[iRow]);
2661	if (value && value < 1.0e30) {
2662	largest = CoinMax(largest, value);
2663	smallest = CoinMin(smallest, value);
2664	}
2665	value = fabs(model2->rowUpper_[iRow]);
2666	if (value && value < 1.0e30) {
2667	largest = CoinMax(largest, value);
2668	smallest = CoinMin(smallest, value);
2669	}
2670	}
2671	double *saveLower = NULL;
2672	double *saveUpper = NULL;
2673	if (largest < -2.01 * smallest) { // switch off for now
2674	// perturb - so switch off standard
2675	model2->setPerturbation(100);
2676	saveLower = new double[numberRows];
2677	CoinMemcpyN(model2->rowLower_, numberRows, saveLower);
2678	saveUpper = new double[numberRows];
2679	CoinMemcpyN(model2->rowUpper_, numberRows, saveUpper);
2680	double *lower = model2->rowLower();
2681	double *upper = model2->rowUpper();
2682	for (iRow = 0; iRow < numberRows; iRow++) {
2683	double lowerValue = lower[iRow], upperValue = upper[iRow];
2684	double value = randomNumberGenerator_.randomDouble();
2685	if (upperValue > lowerValue + primalTolerance_) {
2686	if (lowerValue > -1.0e20 && lowerValue)
2687	lowerValue -= value * 1.0e-4 * fabs(lowerValue);
2688	if (upperValue < 1.0e20 && upperValue)
2689	upperValue += value * 1.0e-4 * fabs(upperValue);
2690	} else if (upperValue > 0.0) {
2691	upperValue -= value * 1.0e-4 * fabs(lowerValue);
2692	lowerValue -= value * 1.0e-4 * fabs(lowerValue);
2693	} else if (upperValue < 0.0) {
2694	upperValue += value * 1.0e-4 * fabs(lowerValue);
2695	lowerValue += value * 1.0e-4 * fabs(lowerValue);
2696	} else {
2697	}
2698	lower[iRow] = lowerValue;
2699	upper[iRow] = upperValue;
2700	}
2701	}
2702	int i;
2703	// Just do this number of passes in Sprint
2704	if (doSprint > 0)
2705	maxSprintPass = options.getExtraInfo(1);
2706	// but if big use to get ratio
2707	double ratio = 3;
2708	#ifdef CLP_USEFUL_PRINTOUT
2709	debugInt[6] = 3;
2710	debugInt[7] = maxSprintPass;
2711	#endif
2712	if (maxSprintPass > 1000) {
2713	ratio = static_cast< double >(maxSprintPass) * 0.0001;
2714	ratio = CoinMax(ratio, 1.1);
2715	maxSprintPass = maxSprintPass % 1000;
2716	#ifdef COIN_DEVELOP
2717	printf("%d passes wanted with ratio of %g\n", maxSprintPass, ratio);
2718	#endif
2719	}
2720	// Just take this number of columns in small problem
2721	int smallNumberColumns = static_cast< int >(CoinMin(ratio * numberRows, static_cast< double >(numberColumns)));
2722	smallNumberColumns = CoinMax(smallNumberColumns, 3000);
2723	smallNumberColumns = CoinMin(smallNumberColumns, numberColumns);
2724	int saveSmallNumber = smallNumberColumns;
2725	bool emergencyMode = false;
2726	//int smallNumberColumns = CoinMin(12*numberRows/10,numberColumns);
2727	//smallNumberColumns = CoinMax(smallNumberColumns,3000);
2728	//smallNumberColumns = CoinMax(smallNumberColumns,numberRows+1000);
2729	// redo as may have changed
2730	columnLower = model2->columnLower();
2731	columnUpper = model2->columnUpper();
2732	columnSolution = model2->primalColumnSolution();
2733	// Set up initial list
2734	numberSort = 0;
2735	if (numberArtificials) {
2736	numberSort = numberArtificials;
2737	for (i = 0; i < numberSort; i++)
2738	sort[i] = i + originalNumberColumns;
2739	}
2740	// put in free
2741	// maybe a solution there already
2742	for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
2743	if (model2->getColumnStatus(iColumn) == basic \|\| (columnLower[iColumn] < -1.0e30 && columnUpper[iColumn] > 1.0e30))
2744	sort[numberSort++] = iColumn;
2745	}
2746	for (iColumn = 0; iColumn < originalNumberColumns; iColumn++) {
2747	if (model2->getColumnStatus(iColumn) != basic) {
2748	if (columnSolution[iColumn] > columnLower[iColumn] && columnSolution[iColumn] < columnUpper[iColumn] && columnSolution[iColumn])
2749	sort[numberSort++] = iColumn;
2750	}
2751	}
2752	numberSort = CoinMin(numberSort, smallNumberColumns);
2753
2754	int numberColumns = model2->numberColumns();
2755	double *fullSolution = model2->primalColumnSolution();
2756
2757	int iPass;
2758	double lastObjective[] = { 1.0e31, 1.0e31 };
2759	// It will be safe to allow dense
2760	model2->setInitialDenseFactorization(true);
2761
2762	// We will be using all rows
2763	int *whichRows = new int[numberRows];
2764	for (iRow = 0; iRow < numberRows; iRow++)
2765	whichRows[iRow] = iRow;
2766	double originalOffset;
2767	model2->getDblParam(ClpObjOffset, originalOffset);
2768	int totalIterations = 0;
2769	double lastSumArtificials = COIN_DBL_MAX;
2770	int originalMaxSprintPass = maxSprintPass;
2771	maxSprintPass = 20; // so we do that many if infeasible
2772	for (iPass = 0; iPass < maxSprintPass; iPass++) {
2773	//printf("Bug until submodel new version\n");
2774	//CoinSort_2(sort,sort+numberSort,weight);
2775	// Create small problem
2776	ClpSimplex small(model2, numberRows, whichRows, numberSort, sort);
2777	small.setPerturbation(model2->perturbation());
2778	small.setInfeasibilityCost(model2->infeasibilityCost());
2779	if (model2->factorizationFrequency() == 200) {
2780	// User did not touch preset
2781	small.defaultFactorizationFrequency();
2782	}
2783	// now see what variables left out do to row solution
2784	double *rowSolution = model2->primalRowSolution();
2785	double *sumFixed = new double[numberRows];
2786	CoinZeroN(sumFixed, numberRows);
2787	int iRow, iColumn;
2788	// zero out ones in small problem
2789	for (iColumn = 0; iColumn < numberSort; iColumn++) {
2790	int kColumn = sort[iColumn];
2791	fullSolution[kColumn] = 0.0;
2792	}
2793	// Get objective offset
2794	const double *objective = model2->objective();
2795	double offset = 0.0;
2796	for (iColumn = 0; iColumn < originalNumberColumns; iColumn++)
2797	offset += fullSolution[iColumn] * objective[iColumn];
2798	#if 0
2799	// Set artificials to zero if first time close to zero
2800	for (iColumn = originalNumberColumns; iColumn < numberColumns; iColumn++) {
2801	if (fullSolution[iColumn]<primalTolerance_&&objective[iColumn]==penalty) {
2802	model2->objective()[iColumn]=2.0*penalty;
2803	fullSolution[iColumn]=0.0;
2804	}
2805	}
2806	#endif
2807	small.setDblParam(ClpObjOffset, originalOffset - offset);
2808	model2->clpMatrix()->times(1.0, fullSolution, sumFixed);
2809
2810	double *lower = small.rowLower();
2811	double *upper = small.rowUpper();
2812	for (iRow = 0; iRow < numberRows; iRow++) {
2813	if (lower[iRow] > -1.0e50)
2814	lower[iRow] -= sumFixed[iRow];
2815	if (upper[iRow] < 1.0e50)
2816	upper[iRow] -= sumFixed[iRow];
2817	rowSolution[iRow] -= sumFixed[iRow];
2818	}
2819	delete[] sumFixed;
2820	// Solve
2821	if (interrupt)
2822	currentModel = &small;
2823	small.defaultFactorizationFrequency();
2824	if (emergencyMode) {
2825	// not much happening so big model
2826	int options = small.moreSpecialOptions();
2827	small.setMoreSpecialOptions(options \| 1048576);
2828	small.setMaximumIterations(1000400);
2829	small.setPerturbation(100);
2830	}
2831	if (dynamic_cast< ClpPackedMatrix * >(matrix_)) {
2832	// See if original wanted vector
2833	ClpPackedMatrix clpMatrixO = dynamic_cast< ClpPackedMatrix >(matrix_);
2834	ClpMatrixBase *matrix = small.clpMatrix();
2835	if (dynamic_cast< ClpPackedMatrix * >(matrix) && clpMatrixO->wantsSpecialColumnCopy()) {
2836	ClpPackedMatrix clpMatrix = dynamic_cast< ClpPackedMatrix >(matrix);
2837	clpMatrix->makeSpecialColumnCopy();
2838	small.primal(1);
2839	clpMatrix->releaseSpecialColumnCopy();
2840	} else {
2841	#if 1
2842	#ifdef ABC_INHERIT
2843	//small.writeMps("try.mps");
2844	if (iPass \|\| !numberArtificials)
2845	small.dealWithAbc(1, 1);
2846	else
2847	small.dealWithAbc(0, 0);
2848	#else
2849	if (iPass \|\| !numberArtificials)
2850	small.primal(1);
2851	else
2852	small.dual(0);
2853	#endif
2854	if (emergencyMode) {
2855	if (small.problemStatus() == 3)
2856	small.setProblemStatus(0);
2857	smallNumberColumns = saveSmallNumber;
2858	emergencyMode = false;
2859	double *temp = new double[numberRows];
2860	memset(temp, 0, numberRows * sizeof(double));
2861	double *solution = small.primalColumnSolution();
2862	small.matrix()->times(solution, temp);
2863	double sumInf = 0.0;
2864	double *lower = small.rowLower();
2865	double *upper = small.rowUpper();
2866	for (int iRow = 0; iRow < numberRows; iRow++) {
2867	if (temp[iRow] > upper[iRow])
2868	sumInf += temp[iRow] - upper[iRow];
2869	else if (temp[iRow] < lower[iRow])
2870	sumInf += lower[iRow] - temp[iRow];
2871	}
2872	printf("row inf %g\n", sumInf);
2873	sumInf = 0.0;
2874	lower = small.columnLower();
2875	upper = small.columnUpper();
2876	for (int iColumn = 0; iColumn < small.numberColumns(); iColumn++) {
2877	if (solution[iColumn] > upper[iColumn])
2878	sumInf += solution[iColumn] - upper[iColumn];
2879	else if (solution[iColumn] < lower[iColumn])
2880	sumInf += lower[iColumn] - solution[iColumn];
2881	}
2882	printf("column inf %g\n", sumInf);
2883	delete[] temp;
2884	}
2885	if (small.problemStatus()) {
2886	int numberIterations = small.numberIterations();
2887	small.dual(0);
2888	small.setNumberIterations(small.numberIterations() + numberIterations);
2889	}
2890	#else
2891	int numberColumns = small.numberColumns();
2892	int numberRows = small.numberRows();
2893	// Use dual region
2894	double *rhs = small.dualRowSolution();
2895	int whichRow = new int[3 numberRows];
2896	int whichColumn = new int[2 numberColumns];
2897	int nBound;
2898	ClpSimplex small2 = ((ClpSimplexOther )(&small))->crunch(rhs, whichRow, whichColumn, nBound, false, false);
2899	if (small2) {
2900	#ifdef ABC_INHERIT
2901	small2->dealWithAbc(1, 1);
2902	#else
2903	small.primal(1);
2904	#endif
2905	if (small2->problemStatus() == 0) {
2906	small.setProblemStatus(0);
2907	((ClpSimplexOther )(&small))->afterCrunch(small2, whichRow, whichColumn, nBound);
2908	} else {
2909	#ifdef ABC_INHERIT
2910	small2->dealWithAbc(1, 1);
2911	#else
2912	small.primal(1);
2913	#endif
2914	if (small2->problemStatus())
2915	small.primal(1);
2916	}
2917	delete small2;
2918	} else {
2919	small.primal(1);
2920	}
2921	delete[] whichRow;
2922	delete[] whichColumn;
2923	#endif
2924	}
2925	} else {
2926	small.primal(1);
2927	}
2928	int smallIterations = small.numberIterations();
2929	totalIterations += smallIterations;
2930	if (2 * smallIterations < CoinMin(numberRows, 1000) && iPass) {
2931	int oldNumber = smallNumberColumns;
2932	if (smallIterations < 100)
2933	smallNumberColumns *= 1.2;
2934	else
2935	smallNumberColumns *= 1.1;
2936	smallNumberColumns = CoinMin(smallNumberColumns, numberColumns);
2937	if (smallIterations < 200 && smallNumberColumns > 2 * saveSmallNumber) {
2938	// try kicking it
2939	emergencyMode = true;
2940	smallNumberColumns = numberColumns;
2941	}
2942	// smallNumberColumns = CoinMin(smallNumberColumns, 3*saveSmallNumber);
2943	char line[100];
2944	sprintf(line, "sample size increased from %d to %d",
2945	oldNumber, smallNumberColumns);
2946	handler_->message(CLP_GENERAL, messages_)
2947	<< line
2948	<< CoinMessageEol;
2949	}
2950	// move solution back
2951	const double *solution = small.primalColumnSolution();
2952	for (iColumn = 0; iColumn < numberSort; iColumn++) {
2953	int kColumn = sort[iColumn];
2954	model2->setColumnStatus(kColumn, small.getColumnStatus(iColumn));
2955	fullSolution[kColumn] = solution[iColumn];
2956	}
2957	for (iRow = 0; iRow < numberRows; iRow++)
2958	model2->setRowStatus(iRow, small.getRowStatus(iRow));
2959	CoinMemcpyN(small.primalRowSolution(),
2960	numberRows, model2->primalRowSolution());
2961	double sumArtificials = 0.0;
2962	for (i = 0; i < numberArtificials; i++)
2963	sumArtificials += fullSolution[i + originalNumberColumns];
2964	if (sumArtificials && iPass > 5 && sumArtificials >= lastSumArtificials) {
2965	// increase costs
2966	double *cost = model2->objective() + originalNumberColumns;
2967	double newCost = CoinMin(1.0e10, cost[0] * 1.5);
2968	for (i = 0; i < numberArtificials; i++)
2969	cost[i] = newCost;
2970	}
2971	lastSumArtificials = sumArtificials;
2972	// get reduced cost for large problem
2973	double *djs = model2->dualColumnSolution();
2974	CoinMemcpyN(model2->objective(), numberColumns, djs);
2975	model2->clpMatrix()->transposeTimes(-1.0, small.dualRowSolution(), djs);
2976	int numberNegative = 0;
2977	double sumNegative = 0.0;
2978	// now massage weight so all basic in plus good djs
2979	// first count and do basic
2980	numberSort = 0;
2981	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
2982	double dj = djs[iColumn] * optimizationDirection_;
2983	double value = fullSolution[iColumn];
2984	if (model2->getColumnStatus(iColumn) == ClpSimplex::basic) {
2985	sort[numberSort++] = iColumn;
2986	} else if (dj < -dualTolerance_ && value < columnUpper[iColumn]) {
2987	numberNegative++;
2988	sumNegative -= dj;
2989	} else if (dj > dualTolerance_ && value > columnLower[iColumn]) {
2990	numberNegative++;
2991	sumNegative += dj;
2992	}
2993	}
2994	handler_->message(CLP_SPRINT, messages_)
2995	<< iPass + 1 << small.numberIterations() << small.objectiveValue() << sumNegative
2996	<< numberNegative
2997	<< CoinMessageEol;
2998	if (sumArtificials < 1.0e-8 && originalMaxSprintPass >= 0) {
2999	maxSprintPass = iPass + originalMaxSprintPass;
3000	originalMaxSprintPass = -1;
3001	}
3002	if (iPass > 20)
3003	sumArtificials = 0.0;
3004	if ((small.objectiveValue() * optimizationDirection_ > lastObjective[1] - 1.0e-7 && iPass > 15 && sumArtificials < 1.0e-8 && maxSprintPass < 200) \|\| (!small.numberIterations() && iPass) \|\| iPass == maxSprintPass - 1 \|\| small.status() == 3) {
3005
3006	break; // finished
3007	} else {
3008	lastObjective[1] = lastObjective[0];
3009	lastObjective[0] = small.objectiveValue() * optimizationDirection_;
3010	double tolerance;
3011	double averageNegDj = sumNegative / static_cast< double >(numberNegative + 1);
3012	if (numberNegative + numberSort > smallNumberColumns && false)
3013	tolerance = -dualTolerance_;
3014	else
3015	tolerance = 10.0 * averageNegDj;
3016	if (emergencyMode)
3017	tolerance = 1.0e100;
3018	int saveN = numberSort;
3019	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
3020	double dj = djs[iColumn] * optimizationDirection_;
3021	double value = fullSolution[iColumn];
3022	if (model2->getColumnStatus(iColumn) != ClpSimplex::basic) {
3023	if (dj < -dualTolerance_ && value < columnUpper[iColumn])
3024	dj = dj;
3025	else if (dj > dualTolerance_ && value > columnLower[iColumn])
3026	dj = -dj;
3027	else if (columnUpper[iColumn] > columnLower[iColumn])
3028	dj = fabs(dj);
3029	else
3030	dj = 1.0e50;
3031	if (dj < tolerance) {
3032	weight[numberSort] = dj;
3033	sort[numberSort++] = iColumn;
3034	}
3035	}
3036	}
3037	// sort
3038	CoinSort_2(weight + saveN, weight + numberSort, sort + saveN);
3039	if (numberSort < smallNumberColumns)
3040	printf("using %d columns not %d\n", numberSort, smallNumberColumns);
3041	numberSort = CoinMin(smallNumberColumns, numberSort);
3042	// try singletons
3043	char *markX = new char[numberColumns];
3044	memset(markX, 0, numberColumns);
3045	for (int i = 0; i < numberSort; i++)
3046	markX[sort[i]] = 1;
3047	int n = numberSort;
3048	for (int i = 0; i < numberColumns; i++) {
3049	if (columnLength[i] == 1 && !markX[i])
3050	sort[numberSort++] = i;
3051	}
3052	if (n < numberSort)
3053	printf("%d slacks added\n", numberSort - n);
3054	delete[] markX;
3055	}
3056	}
3057	if (interrupt)
3058	currentModel = model2;
3059	for (i = 0; i < numberArtificials; i++)
3060	sort[i] = i + originalNumberColumns;
3061	model2->deleteColumns(numberArtificials, sort);
3062	if (network) {
3063	int iRow = numberRows - 1;
3064	model2->deleteRows(1, &iRow);
3065	}
3066	delete[] weight;
3067	delete[] sort;
3068	delete[] whichRows;
3069	if (saveLower) {
3070	// unperturb and clean
3071	for (iRow = 0; iRow < numberRows; iRow++) {
3072	model2->rowLower_[iRow] = saveLower[iRow];
3073	model2->rowUpper_[iRow] = saveUpper[iRow];
3074	}
3075	delete[] saveLower;
3076	delete[] saveUpper;
3077	}
3078	#ifdef ABC_INHERIT
3079	model2->dealWithAbc(1, 1);
3080	#else
3081	model2->primal(1);
3082	#endif
3083	model2->setPerturbation(savePerturbation);
3084	if (model2 != originalModel2) {
3085	originalModel2->moveInfo(*model2);
3086	delete model2;
3087	model2 = originalModel2;
3088	}
3089	time2 = CoinCpuTime();
3090	timeCore = time2 - timeX;
3091	handler_->message(CLP_INTERVAL_TIMING, messages_)
3092	<< "Sprint" << timeCore << time2 - time1
3093	<< CoinMessageEol;
3094	timeX = time2;
3095	model2->setNumberIterations(model2->numberIterations() + totalIterations);
3096	} else if (method == ClpSolve::useBarrier \|\| method == ClpSolve::useBarrierNoCross) {
3097	if (presolve == ClpSolve::presolveOn) {
3098	int numberInfeasibilities = model2->tightenPrimalBounds(0.0, 0);
3099	if (numberInfeasibilities) {
3100	handler_->message(CLP_INFEASIBLE, messages_)
3101	<< CoinMessageEol;
3102	delete model2;
3103	model2 = this;
3104	presolve = ClpSolve::presolveOff;
3105	}
3106	}
3107	#ifndef SLIM_CLP
3108	//printf("***** experimental pretty crude barrier\n");
3109	//#define SAVEIT 2
3110	#ifndef SAVEIT
3111	#define BORROW
3112	#endif
3113	#ifdef BORROW
3114	ClpInterior barrier;
3115	barrier.borrowModel(*model2);
3116	#else
3117	ClpInterior barrier(*model2);
3118	#endif
3119	if (interrupt)
3120	currentModel2 = &barrier;
3121	if (barrier.numberRows() + barrier.numberColumns() > 10000)
3122	barrier.setMaximumBarrierIterations(1000);
3123	int barrierOptions = options.getSpecialOption(4);
3124	int aggressiveGamma = 0;
3125	bool presolveInCrossover = false;
3126	bool scale = false;
3127	bool doKKT = false;
3128	bool forceFixing = false;
3129	int speed = 0;
3130	if (barrierOptions & 16) {
3131	barrierOptions &= ~16;
3132	doKKT = true;
3133	}
3134	if (barrierOptions & (32 + 64 + 128)) {
3135	aggressiveGamma = (barrierOptions & (32 + 64 + 128)) >> 5;
3136	barrierOptions &= ~(32 + 64 + 128);
3137	}
3138	if (barrierOptions & 256) {
3139	barrierOptions &= ~256;
3140	presolveInCrossover = true;
3141	}
3142	if (barrierOptions & 512) {
3143	barrierOptions &= ~512;
3144	forceFixing = true;
3145	}
3146	if (barrierOptions & 1024) {
3147	barrierOptions &= ~1024;
3148	barrier.setProjectionTolerance(1.0e-9);
3149	}
3150	if (barrierOptions & (2048 \| 4096)) {
3151	speed = (barrierOptions & (2048 \| 4096)) >> 11;
3152	barrierOptions &= ~(2048 \| 4096);
3153	}
3154	if (barrierOptions & 8) {
3155	barrierOptions &= ~8;
3156	scale = true;
3157	}
3158	// If quadratic force KKT
3159	if (quadraticObj) {
3160	doKKT = true;
3161	}
3162	switch (barrierOptions) {
3163	case 0:
3164	default:
3165	if (!doKKT) {
3166	ClpCholeskyBase *cholesky = new ClpCholeskyBase(options.getExtraInfo(1));
3167	cholesky->setIntegerParameter(0, speed);
3168	barrier.setCholesky(cholesky);
3169	} else {
3170	ClpCholeskyBase *cholesky = new ClpCholeskyBase();
3171	cholesky->setKKT(true);
3172	barrier.setCholesky(cholesky);
3173	}
3174	break;
3175	case 1:
3176	if (!doKKT) {
3177	ClpCholeskyDense *cholesky = new ClpCholeskyDense();
3178	barrier.setCholesky(cholesky);
3179	} else {
3180	ClpCholeskyDense *cholesky = new ClpCholeskyDense();
3181	cholesky->setKKT(true);
3182	barrier.setCholesky(cholesky);
3183	}
3184	break;
3185	#ifdef COIN_HAS_WSMP
3186	case 2: {
3187	if (!doKKT) {
3188	ClpCholeskyWssmp *cholesky = new ClpCholeskyWssmp(CoinMax(100, model2->numberRows() / 10));
3189	barrier.setCholesky(cholesky);
3190	} else {
3191	ClpCholeskyWssmpKKT *cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
3192	barrier.setCholesky(cholesky);
3193	}
3194	} break;
3195	case 3:
3196	if (!doKKT) {
3197	ClpCholeskyWssmp *cholesky = new ClpCholeskyWssmp();
3198	barrier.setCholesky(cholesky);
3199	} else {
3200	ClpCholeskyWssmpKKT *cholesky = new ClpCholeskyWssmpKKT(CoinMax(100, model2->numberRows() / 10));
3201	barrier.setCholesky(cholesky);
3202	}
3203	break;
3204	#endif
3205	#ifdef UFL_BARRIER
3206	case 4:
3207	if (!doKKT) {
3208	ClpCholeskyUfl *cholesky = new ClpCholeskyUfl(options.getExtraInfo(1));
3209	barrier.setCholesky(cholesky);
3210	} else {
3211	ClpCholeskyUfl *cholesky = new ClpCholeskyUfl();
3212	cholesky->setKKT(true);
3213	barrier.setCholesky(cholesky);
3214	}
3215	break;
3216	#endif
3217	#ifdef TAUCS_BARRIER
3218	case 5: {
3219	ClpCholeskyTaucs *cholesky = new ClpCholeskyTaucs();
3220	barrier.setCholesky(cholesky);
3221	assert(!doKKT);
3222	} break;
3223	#endif
3224	#ifdef COIN_HAS_MUMPS
3225	case 6: {
3226	if (!doKKT) {
3227	int logLevel = 0;
3228	if (barrier.logLevel() > 3) {
3229	logLevel = (barrier.logLevel() > 4) ? 2 : 1;
3230	}
3231	ClpCholeskyMumps *cholesky = new ClpCholeskyMumps(-1, logLevel);
3232	barrier.setCholesky(cholesky);
3233	} else {
3234	ClpCholeskyMumps *cholesky = new ClpCholeskyMumps();
3235	cholesky->setKKT(true);
3236	barrier.setCholesky(cholesky);
3237	}
3238	} break;
3239	#endif
3240	#ifdef PARDISO_BARRIER
3241	case 7: {
3242	ClpCholeskyPardiso *cholesky = new ClpCholeskyPardiso();
3243	barrier.setCholesky(cholesky);
3244	assert(!doKKT);
3245	} break;
3246	#endif
3247	}
3248	int numberRows = model2->numberRows();
3249	int numberColumns = model2->numberColumns();
3250	int saveMaxIts = model2->maximumIterations();
3251	if (saveMaxIts < 1000) {
3252	barrier.setMaximumBarrierIterations(saveMaxIts);
3253	model2->setMaximumIterations(10000000);
3254	}
3255	#ifndef SAVEIT
3256	//barrier.setDiagonalPerturbation(1.0e-25);
3257	if (aggressiveGamma) {
3258	switch (aggressiveGamma) {
3259	case 1:
3260	barrier.setGamma(1.0e-5);
3261	barrier.setDelta(1.0e-5);
3262	break;
3263	case 2:
3264	barrier.setGamma(1.0e-7);
3265	break;
3266	case 3:
3267	barrier.setDelta(1.0e-5);
3268	break;
3269	case 4:
3270	barrier.setGamma(1.0e-3);
3271	barrier.setDelta(1.0e-3);
3272	break;
3273	case 5:
3274	barrier.setGamma(1.0e-3);
3275	break;
3276	case 6:
3277	barrier.setDelta(1.0e-3);
3278	break;
3279	}
3280	}
3281	if (scale)
3282	barrier.scaling(1);
3283	else
3284	barrier.scaling(0);
3285	barrier.primalDual();
3286	#elif SAVEIT == 1
3287	barrier.primalDual();
3288	#else
3289	model2->restoreModel("xx.save");
3290	// move solutions
3291	CoinMemcpyN(model2->primalRowSolution(),
3292	numberRows, barrier.primalRowSolution());
3293	CoinMemcpyN(model2->dualRowSolution(),
3294	numberRows, barrier.dualRowSolution());
3295	CoinMemcpyN(model2->primalColumnSolution(),
3296	numberColumns, barrier.primalColumnSolution());
3297	CoinMemcpyN(model2->dualColumnSolution(),
3298	numberColumns, barrier.dualColumnSolution());
3299	#endif
3300	time2 = CoinCpuTime();
3301	timeCore = time2 - timeX;
3302	handler_->message(CLP_INTERVAL_TIMING, messages_)
3303	<< "Barrier" << timeCore << time2 - time1
3304	<< CoinMessageEol;
3305	timeX = time2;
3306	int maxIts = barrier.maximumBarrierIterations();
3307	int barrierStatus = barrier.status();
3308	double gap = barrier.complementarityGap();
3309	// get which variables are fixed
3310	double *saveLower = NULL;
3311	double *saveUpper = NULL;
3312	ClpPresolve pinfo2;
3313	ClpSimplex *saveModel2 = NULL;
3314	bool extraPresolve = false;
3315	int numberFixed = barrier.numberFixed();
3316	if (numberFixed) {
3317	int numberRows = barrier.numberRows();
3318	int numberColumns = barrier.numberColumns();
3319	int numberTotal = numberRows + numberColumns;
3320	saveLower = new double[numberTotal];
3321	saveUpper = new double[numberTotal];
3322	CoinMemcpyN(barrier.columnLower(), numberColumns, saveLower);
3323	CoinMemcpyN(barrier.rowLower(), numberRows, saveLower + numberColumns);
3324	CoinMemcpyN(barrier.columnUpper(), numberColumns, saveUpper);
3325	CoinMemcpyN(barrier.rowUpper(), numberRows, saveUpper + numberColumns);
3326	}
3327	if (((numberFixed * 20 > barrier.numberRows() && numberFixed > 5000) \|\| forceFixing) && presolveInCrossover) {
3328	// may as well do presolve
3329	if (!forceFixing) {
3330	barrier.fixFixed();
3331	} else {
3332	// Fix
3333	int n = barrier.numberColumns();
3334	double *lower = barrier.columnLower();
3335	double *upper = barrier.columnUpper();
3336	double *solution = barrier.primalColumnSolution();
3337	int nFix = 0;
3338	for (int i = 0; i < n; i++) {
3339	if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
3340	double value = solution[i];
3341	if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
3342	solution[i] = lower[i];
3343	upper[i] = lower[i];
3344	nFix++;
3345	} else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
3346	solution[i] = upper[i];
3347	lower[i] = upper[i];
3348	nFix++;
3349	}
3350	}
3351	}
3352	#ifdef CLP_INVESTIGATE
3353	printf("%d columns fixed\n", nFix);
3354	#endif
3355	int nr = barrier.numberRows();
3356	lower = barrier.rowLower();
3357	upper = barrier.rowUpper();
3358	solution = barrier.primalRowSolution();
3359	nFix = 0;
3360	for (int i = 0; i < nr; i++) {
3361	if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
3362	double value = solution[i];
3363	if (value < lower[i] + 1.0e-6 && value - lower[i] < upper[i] - value) {
3364	solution[i] = lower[i];
3365	upper[i] = lower[i];
3366	nFix++;
3367	} else if (value > upper[i] - 1.0e-6 && value - lower[i] > upper[i] - value) {
3368	solution[i] = upper[i];
3369	lower[i] = upper[i];
3370	nFix++;
3371	}
3372	}
3373	}
3374	#ifdef CLP_INVESTIGATE
3375	printf("%d row slacks fixed\n", nFix);
3376	#endif
3377	}
3378	saveModel2 = model2;
3379	extraPresolve = true;
3380	} else if (numberFixed) {
3381	// Set fixed to bounds (may have restored earlier solution)
3382	if (!forceFixing) {
3383	barrier.fixFixed(false);
3384	} else {
3385	// Fix
3386	int n = barrier.numberColumns();
3387	double *lower = barrier.columnLower();
3388	double *upper = barrier.columnUpper();
3389	double *solution = barrier.primalColumnSolution();
3390	int nFix = 0;
3391	for (int i = 0; i < n; i++) {
3392	if (barrier.fixedOrFree(i) && lower[i] < upper[i]) {
3393	double value = solution[i];
3394	if (value < lower[i] + 1.0e-8 && value - lower[i] < upper[i] - value) {
3395	solution[i] = lower[i];
3396	upper[i] = lower[i];
3397	nFix++;
3398	} else if (value > upper[i] - 1.0e-8 && value - lower[i] > upper[i] - value) {
3399	solution[i] = upper[i];
3400	lower[i] = upper[i];
3401	nFix++;
3402	} else {
3403	//printf("fixcol %d %g <= %g <= %g\n",
3404	// i,lower[i],solution[i],upper[i]);
3405	}
3406	}
3407	}
3408	#ifdef CLP_INVESTIGATE
3409	printf("%d columns fixed\n", nFix);
3410	#endif
3411	int nr = barrier.numberRows();
3412	lower = barrier.rowLower();
3413	upper = barrier.rowUpper();
3414	solution = barrier.primalRowSolution();
3415	nFix = 0;
3416	for (int i = 0; i < nr; i++) {
3417	if (barrier.fixedOrFree(i + n) && lower[i] < upper[i]) {
3418	double value = solution[i];
3419	if (value < lower[i] + 1.0e-5 && value - lower[i] < upper[i] - value) {
3420	solution[i] = lower[i];
3421	upper[i] = lower[i];
3422	nFix++;
3423	} else if (value > upper[i] - 1.0e-5 && value - lower[i] > upper[i] - value) {
3424	solution[i] = upper[i];
3425	lower[i] = upper[i];
3426	nFix++;
3427	} else {
3428	//printf("fixrow %d %g <= %g <= %g\n",
3429	// i,lower[i],solution[i],upper[i]);
3430	}
3431	}
3432	}
3433	#ifdef CLP_INVESTIGATE
3434	printf("%d row slacks fixed\n", nFix);
3435	#endif
3436	}
3437	}
3438	#ifdef BORROW
3439	int saveNumberIterations = barrier.numberIterations();
3440	barrier.returnModel(*model2);
3441	double *rowPrimal = new double[numberRows];
3442	double *columnPrimal = new double[numberColumns];
3443	double *rowDual = new double[numberRows];
3444	double *columnDual = new double[numberColumns];
3445	// move solutions other way
3446	CoinMemcpyN(model2->primalRowSolution(),
3447	numberRows, rowPrimal);
3448	CoinMemcpyN(model2->dualRowSolution(),
3449	numberRows, rowDual);
3450	CoinMemcpyN(model2->primalColumnSolution(),
3451	numberColumns, columnPrimal);
3452	CoinMemcpyN(model2->dualColumnSolution(),
3453	numberColumns, columnDual);
3454	#else
3455	double *rowPrimal = barrier.primalRowSolution();
3456	double *columnPrimal = barrier.primalColumnSolution();
3457	double *rowDual = barrier.dualRowSolution();
3458	double *columnDual = barrier.dualColumnSolution();
3459	// move solutions
3460	CoinMemcpyN(rowPrimal,
3461	numberRows, model2->primalRowSolution());
3462	CoinMemcpyN(rowDual,
3463	numberRows, model2->dualRowSolution());
3464	CoinMemcpyN(columnPrimal,
3465	numberColumns, model2->primalColumnSolution());
3466	CoinMemcpyN(columnDual,
3467	numberColumns, model2->dualColumnSolution());
3468	#endif
3469	if (saveModel2) {
3470	// do presolve
3471	model2 = pinfo2.presolvedModel(*model2, dblParam_[ClpPresolveTolerance],
3472	false, 5, true);
3473	if (!model2) {
3474	model2 = saveModel2;
3475	saveModel2 = NULL;
3476	int numberRows = model2->numberRows();
3477	int numberColumns = model2->numberColumns();
3478	CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3479	CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3480	delete[] saveLower;
3481	CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3482	CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3483	delete[] saveUpper;
3484	saveLower = NULL;
3485	saveUpper = NULL;
3486	}
3487	}
3488	if (method == ClpSolve::useBarrier \|\| barrierStatus < 0) {
3489	if (maxIts && barrierStatus < 4 && !quadraticObj) {
3490	//printf("***** crossover - needs more thought on difficult models\n");
3491	#if SAVEIT == 1
3492	model2->ClpSimplex::saveModel("xx.save");
3493	#endif
3494	// make sure no status left
3495	model2->createStatus();
3496	// solve
3497	if (!forceFixing)
3498	model2->setPerturbation(100);
3499	if (model2->factorizationFrequency() == 200) {
3500	// User did not touch preset
3501	model2->defaultFactorizationFrequency();
3502	}
3503	#if 1 //ndef ABC_INHERIT //#if 1 \
3504	// throw some into basis
3505	if (!forceFixing) {
3506	int numberRows = model2->numberRows();
3507	int numberColumns = model2->numberColumns();
3508	double *dsort = new double[numberColumns];
3509	int *sort = new int[numberColumns];
3510	int n = 0;
3511	const double *columnLower = model2->columnLower();
3512	const double *columnUpper = model2->columnUpper();
3513	double *primalSolution = model2->primalColumnSolution();
3514	const double *dualSolution = model2->dualColumnSolution();
3515	double tolerance = 10.0 * primalTolerance_;
3516	int i;
3517	for (i = 0; i < numberRows; i++)
3518	model2->setRowStatus(i, superBasic);
3519	for (i = 0; i < numberColumns; i++) {
3520	double distance = CoinMin(columnUpper[i] - primalSolution[i],
3521	primalSolution[i] - columnLower[i]);
3522	if (distance > tolerance) {
3523	if (fabs(dualSolution[i]) < 1.0e-5)
3524	distance *= 100.0;
3525	dsort[n] = -distance;
3526	sort[n++] = i;
3527	model2->setStatus(i, superBasic);
3528	} else if (distance > primalTolerance_) {
3529	model2->setStatus(i, superBasic);
3530	} else if (primalSolution[i] <= columnLower[i] + primalTolerance_) {
3531	model2->setStatus(i, atLowerBound);
3532	primalSolution[i] = columnLower[i];
3533	} else {
3534	model2->setStatus(i, atUpperBound);
3535	primalSolution[i] = columnUpper[i];
3536	}
3537	}
3538	CoinSort_2(dsort, dsort + n, sort);
3539	n = CoinMin(numberRows, n);
3540	for (i = 0; i < n; i++) {
3541	int iColumn = sort[i];
3542	model2->setStatus(iColumn, basic);
3543	}
3544	delete[] sort;
3545	delete[] dsort;
3546	// model2->allSlackBasis();
3547	if (gap < 1.0e-3 * static_cast< double >(numberRows + numberColumns)) {
3548	if (saveUpper) {
3549	int numberRows = model2->numberRows();
3550	int numberColumns = model2->numberColumns();
3551	CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3552	CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3553	CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3554	CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3555	delete[] saveLower;
3556	delete[] saveUpper;
3557	saveLower = NULL;
3558	saveUpper = NULL;
3559	}
3560	//int numberRows = model2->numberRows();
3561	//int numberColumns = model2->numberColumns();
3562	#ifdef ABC_INHERIT
3563	model2->checkSolution(0);
3564	printf("%d primal infeasibilities summing to %g\n",
3565	model2->numberPrimalInfeasibilities(),
3566	model2->sumPrimalInfeasibilities());
3567	model2->dealWithAbc(1, 1);
3568	}
3569	}
3570	#else
3571	// just primal values pass
3572	double saveScale = model2->objectiveScale();
3573	model2->setObjectiveScale(1.0e-3);
3574	model2->primal(2);
3575	model2->setObjectiveScale(saveScale);
3576	// save primal solution and copy back dual
3577	CoinMemcpyN(model2->primalRowSolution(),
3578	numberRows, rowPrimal);
3579	CoinMemcpyN(rowDual,
3580	numberRows, model2->dualRowSolution());
3581	CoinMemcpyN(model2->primalColumnSolution(),
3582	numberColumns, columnPrimal);
3583	CoinMemcpyN(columnDual,
3584	numberColumns, model2->dualColumnSolution());
3585	//model2->primal(1);
3586	// clean up reduced costs and flag variables
3587	{
3588	double *dj = model2->dualColumnSolution();
3589	double *cost = model2->objective();
3590	double *saveCost = new double[numberColumns];
3591	CoinMemcpyN(cost, numberColumns, saveCost);
3592	double *saveLower = new double[numberColumns];
3593	double *lower = model2->columnLower();
3594	CoinMemcpyN(lower, numberColumns, saveLower);
3595	double *saveUpper = new double[numberColumns];
3596	double *upper = model2->columnUpper();
3597	CoinMemcpyN(upper, numberColumns, saveUpper);
3598	int i;
3599	double tolerance = 10.0 * dualTolerance_;
3600	for (i = 0; i < numberColumns; i++) {
3601	if (model2->getStatus(i) == basic) {
3602	dj[i] = 0.0;
3603	} else if (model2->getStatus(i) == atLowerBound) {
3604	if (optimizationDirection_ * dj[i] < tolerance) {
3605	if (optimizationDirection_ * dj[i] < 0.0) {
3606	//if (dj[i]<-1.0e-3)
3607	//printf("bad dj at lb %d %g\n",i,dj[i]);
3608	cost[i] -= dj[i];
3609	dj[i] = 0.0;
3610	}
3611	} else {
3612	upper[i] = lower[i];
3613	}
3614	} else if (model2->getStatus(i) == atUpperBound) {
3615	if (optimizationDirection_ * dj[i] > tolerance) {
3616	if (optimizationDirection_ * dj[i] > 0.0) {
3617	//if (dj[i]>1.0e-3)
3618	//printf("bad dj at ub %d %g\n",i,dj[i]);
3619	cost[i] -= dj[i];
3620	dj[i] = 0.0;
3621	}
3622	} else {
3623	lower[i] = upper[i];
3624	}
3625	}
3626	}
3627	// just dual values pass
3628	//model2->setLogLevel(63);
3629	//model2->setFactorizationFrequency(1);
3630	if (!gap)
3631	model2->dual(2);
3632	CoinMemcpyN(saveCost, numberColumns, cost);
3633	delete[] saveCost;
3634	CoinMemcpyN(saveLower, numberColumns, lower);
3635	delete[] saveLower;
3636	CoinMemcpyN(saveUpper, numberColumns, upper);
3637	delete[] saveUpper;
3638	}
3639	}
3640	// and finish
3641	// move solutions
3642	CoinMemcpyN(rowPrimal,
3643	numberRows, model2->primalRowSolution());
3644	CoinMemcpyN(columnPrimal,
3645	numberColumns, model2->primalColumnSolution());
3646	}
3647	double saveScale = model2->objectiveScale();
3648	model2->setObjectiveScale(1.0e-3);
3649	model2->primal(2);
3650	model2->setObjectiveScale(saveScale);
3651	model2->primal(1);
3652	#endif
3653	#else
3654	// just primal
3655	#ifdef ABC_INHERIT
3656	model2->checkSolution(0);
3657	printf("%d primal infeasibilities summing to %g\n",
3658	model2->numberPrimalInfeasibilities(),
3659	model2->sumPrimalInfeasibilities());
3660	model2->dealWithAbc(1, 1);
3661	#else
3662	model2->primal(1);
3663	#endif
3664	//model2->primal(1);
3665	#endif
3666	} else if (barrierStatus == 4) {
3667	// memory problems
3668	model2->setPerturbation(savePerturbation);
3669	model2->createStatus();
3670	model2->dual();
3671	} else if (maxIts && quadraticObj) {
3672	// make sure no status left
3673	model2->createStatus();
3674	// solve
3675	model2->setPerturbation(100);
3676	model2->reducedGradient(1);
3677	}
3678	}
3679
3680	//model2->setMaximumIterations(saveMaxIts);
3681	#ifdef BORROW
3682	model2->setNumberIterations(model2->numberIterations() + saveNumberIterations);
3683	delete[] rowPrimal;
3684	delete[] columnPrimal;
3685	delete[] rowDual;
3686	delete[] columnDual;
3687	#endif
3688	if (extraPresolve) {
3689	pinfo2.postsolve(true);
3690	delete model2;
3691	model2 = saveModel2;
3692	}
3693	if (saveUpper) {
3694	if (!forceFixing) {
3695	int numberRows = model2->numberRows();
3696	int numberColumns = model2->numberColumns();
3697	CoinMemcpyN(saveLower, numberColumns, model2->columnLower());
3698	CoinMemcpyN(saveLower + numberColumns, numberRows, model2->rowLower());
3699	CoinMemcpyN(saveUpper, numberColumns, model2->columnUpper());
3700	CoinMemcpyN(saveUpper + numberColumns, numberRows, model2->rowUpper());
3701	}
3702	delete[] saveLower;
3703	delete[] saveUpper;
3704	saveLower = NULL;
3705	saveUpper = NULL;
3706	if (method != ClpSolve::useBarrierNoCross)
3707	model2->primal(1);
3708	}
3709	model2->setPerturbation(savePerturbation);
3710	time2 = CoinCpuTime();
3711	timeCore = time2 - timeX;
3712	handler_->message(CLP_INTERVAL_TIMING, messages_)
3713	<< "Crossover" << timeCore << time2 - time1
3714	<< CoinMessageEol;
3715	timeX = time2;
3716	#else
3717	abort();
3718	#endif
3719	} else if (method == ClpSolve::notImplemented) {
3720	printf("done decomposition\n");
3721	} else {
3722	assert(method != ClpSolve::automatic); // later
3723	time2 = 0.0;
3724	}
3725	if (saveMatrix) {
3726	if (model2 == this) {
3727	// delete and replace
3728	delete model2->clpMatrix();
3729	model2->replaceMatrix(saveMatrix);
3730	} else {
3731	delete saveMatrix;
3732	}
3733	}
3734	numberIterations = model2->numberIterations();
3735	finalStatus = model2->status();
3736	int finalSecondaryStatus = model2->secondaryStatus();
3737	if (presolve == ClpSolve::presolveOn) {
3738	int saveLevel = logLevel();
3739	if ((specialOptions_ & 1024) == 0)
3740	setLogLevel(CoinMin(1, saveLevel));
3741	else
3742	setLogLevel(CoinMin(0, saveLevel));
3743	pinfo->postsolve(true);
3744	numberIterations_ = 0;
3745	delete pinfo;
3746	pinfo = NULL;
3747	factorization_->areaFactor(model2->factorization()->adjustedAreaFactor());
3748	time2 = CoinCpuTime();
3749	timePresolve += time2 - timeX;
3750	handler_->message(CLP_INTERVAL_TIMING, messages_)
3751	<< "Postsolve" << time2 - timeX << time2 - time1
3752	<< CoinMessageEol;
3753	timeX = time2;
3754	if (!presolveToFile) {
3755	#if 1 //ndef ABC_INHERIT
3756	delete model2;
3757	#else
3758	if (model2->abcSimplex())
3759	delete model2->abcSimplex();
3760	else
3761	delete model2;
3762	#endif
3763	}
3764	if (interrupt)
3765	currentModel = this;
3766	// checkSolution(); already done by postSolve
3767	setLogLevel(saveLevel);
3768	int oldStatus = problemStatus_;
3769	setProblemStatus(finalStatus);
3770	setSecondaryStatus(finalSecondaryStatus);
3771	/* Code modified so rcode -1 as normal, >0 clean up, 0 say optimal */
3772	int rcode = eventHandler()->event(ClpEventHandler::presolveAfterFirstSolve);
3773	//#define TREAT_AS_OPTIMAL_TOLERANCE 1.0e-4
3774	#ifdef TREAT_AS_OPTIMAL_TOLERANCE
3775	if (rcode == -1 && sumPrimalInfeasibilities_ < TREAT_AS_OPTIMAL_TOLERANCE && sumDualInfeasibilities_ < TREAT_AS_OPTIMAL_TOLERANCE)
3776	rcode = 0;
3777	#endif
3778	if (finalStatus != 3 && rcode < 0 && (finalStatus \|\| oldStatus == -1)) {
3779	double sumPrimal = sumPrimalInfeasibilities_;
3780	double sumDual = sumDualInfeasibilities_;
3781	if (sumDual > 1.0e-6 && sumPrimal > 1.0e-6)
3782	moreSpecialOptions_ &= ~2; // be safe and do final solve
3783	// ignore some parts of solution
3784	if (finalStatus == 1) {
3785	// infeasible
3786	sumDual = 0.0;
3787	} else if (finalStatus == 2) {
3788	sumPrimal = 0.0;
3789	}
3790	int savePerturbation = perturbation();
3791	if (savePerturbation == 50)
3792	setPerturbation(51); // small
3793	if (!finalStatus \|\| finalStatus == 2 \|\| (moreSpecialOptions_ & 2) == 0 \|\| fabs(sumDual) + fabs(sumPrimal) < 1.0e-3) {
3794	if (finalStatus == 2) {
3795	if (sumDual > 1.0e-4) {
3796	// unbounded - get feasible first
3797	double save = optimizationDirection_;
3798	optimizationDirection_ = 0.0;
3799	primal(1);
3800	optimizationDirection_ = save;
3801	}
3802	primal(1);
3803	} else if (finalStatus == 1) {
3804	dual();
3805	} else {
3806	if ((moreSpecialOptions_ & 65536) == 0) {
3807	if (numberRows_ < 10000 \|\| true)
3808	setPerturbation(100); // probably better to perturb after n its
3809	else if (savePerturbation < 100)
3810	setPerturbation(51); // probably better to perturb after n its
3811	}
3812	#ifndef ABC_INHERIT
3813	// use method thought suitable
3814	int numberSuperBasic = 0;
3815	for (int i = 0; i < numberColumns_; i++) {
3816	if (getColumnStatus(i) == superBasic)
3817	numberSuperBasic++;
3818	}
3819	if (sumDual > 1000.0 * sumPrimal \|\| numberSuperBasic) {
3820	primal(1);
3821	} else if (sumPrimal > 1000.0 * sumDual) {
3822	dual();
3823	} else {
3824	if (method != ClpSolve::useDual)
3825	primal(1);
3826	else
3827	dual();
3828	}
3829	#else
3830	dealWithAbc(1, 2, interrupt);
3831	#endif
3832	}
3833	} else {
3834	// just set status
3835	problemStatus_ = finalStatus;
3836	}
3837	setPerturbation(savePerturbation);
3838	numberIterations += numberIterations_;
3839	numberIterations_ = numberIterations;
3840	finalStatus = status();
3841	time2 = CoinCpuTime();
3842	handler_->message(CLP_INTERVAL_TIMING, messages_)
3843	<< "Cleanup" << time2 - timeX << time2 - time1
3844	<< CoinMessageEol;
3845	timeX = time2;
3846	} else if (rcode > 0) { // was >= 0
3847	#ifdef ABC_INHERIT
3848	dealWithAbc(1, 2, true);
3849	#else
3850	primal(1);
3851	#endif
3852	} else {
3853	secondaryStatus_ = finalSecondaryStatus;
3854	}
3855	} else if (model2 != this) {
3856	// not presolved - but different model used (sprint probably)
3857	CoinMemcpyN(model2->primalRowSolution(),
3858	numberRows_, this->primalRowSolution());
3859	CoinMemcpyN(model2->dualRowSolution(),
3860	numberRows_, this->dualRowSolution());
3861	CoinMemcpyN(model2->primalColumnSolution(),
3862	numberColumns_, this->primalColumnSolution());
3863	CoinMemcpyN(model2->dualColumnSolution(),
3864	numberColumns_, this->dualColumnSolution());
3865	CoinMemcpyN(model2->statusArray(),
3866	numberColumns_ + numberRows_, this->statusArray());
3867	objectiveValue_ = model2->objectiveValue_;
3868	numberIterations_ = model2->numberIterations_;
3869	problemStatus_ = model2->problemStatus_;
3870	secondaryStatus_ = model2->secondaryStatus_;
3871	delete model2;
3872	}
3873	if (method != ClpSolve::useBarrierNoCross && method != ClpSolve::useBarrier)
3874	setMaximumIterations(saveMaxIterations);
3875	std::string statusMessage[] = { "Unknown", "Optimal", "PrimalInfeasible", "DualInfeasible", "Stopped",
3876	"Errors", "User stopped" };
3877	assert(finalStatus >= -1 && finalStatus <= 5);
3878	numberIterations_ = numberIterations;
3879	handler_->message(CLP_TIMING, messages_)
3880	<< statusMessage[finalStatus + 1] << objectiveValue() << numberIterations << time2 - time1;
3881	handler_->printing(presolve == ClpSolve::presolveOn)
3882	<< timePresolve;
3883	handler_->printing(timeIdiot != 0.0)
3884	<< timeIdiot;
3885	handler_->message() << CoinMessageEol;
3886	if (interrupt)
3887	signal(SIGINT, saveSignal);
3888	perturbation_ = savePerturbation;
3889	scalingFlag_ = saveScaling;
3890	// If faking objective - put back correct one
3891	if (savedObjective) {
3892	delete objective_;
3893	objective_ = savedObjective;
3894	}
3895	if (options.getSpecialOption(1) == 2 && options.getExtraInfo(1) > 1000000) {
3896	ClpObjective *savedObjective = objective_;
3897	// make up zero objective
3898	double *obj = new double[numberColumns_];
3899	for (int i = 0; i < numberColumns_; i++)
3900	obj[i] = 0.0;
3901	objective_ = new ClpLinearObjective(obj, numberColumns_);
3902	delete[] obj;
3903	primal(1);
3904	delete objective_;
3905	objective_ = savedObjective;
3906	finalStatus = status();
3907	}
3908	eventHandler()->event(ClpEventHandler::presolveEnd);
3909	delete pinfo;
3910	moreSpecialOptions_ = saveMoreOptions;
3911	#ifdef CLP_USEFUL_PRINTOUT
3912	debugInt[23] = numberIterations_;
3913	#endif
3914	return finalStatus;
3915	}
3916	// General solve
3917	int ClpSimplex::initialSolve()
3918	{
3919	// Default so use dual
3920	ClpSolve options;
3921	return initialSolve(options);
3922	}
3923	// General dual solve
3924	int ClpSimplex::initialDualSolve()
3925	{
3926	ClpSolve options;
3927	// Use dual
3928	options.setSolveType(ClpSolve::useDual);
3929	return initialSolve(options);
3930	}
3931	// General primal solve
3932	int ClpSimplex::initialPrimalSolve()
3933	{
3934	ClpSolve options;
3935	// Use primal
3936	options.setSolveType(ClpSolve::usePrimal);
3937	return initialSolve(options);
3938	}
3939	// barrier solve, not to be followed by crossover
3940	int ClpSimplex::initialBarrierNoCrossSolve()
3941	{
3942	ClpSolve options;
3943	// Use primal
3944	options.setSolveType(ClpSolve::useBarrierNoCross);
3945	return initialSolve(options);
3946	}
3947
3948	// General barrier solve
3949	int ClpSimplex::initialBarrierSolve()
3950	{
3951	ClpSolve options;
3952	// Use primal
3953	options.setSolveType(ClpSolve::useBarrier);
3954	return initialSolve(options);
3955	}
3956
3957	// Default constructor
3958	ClpSolve::ClpSolve()
3959	{
3960	method_ = automatic;
3961	presolveType_ = presolveOn;
3962	numberPasses_ = 5;
3963	int i;
3964	for (i = 0; i < 7; i++)
3965	options_[i] = 0;
3966	// say no +-1 matrix
3967	options_[3] = 1;
3968	for (i = 0; i < 7; i++)
3969	extraInfo_[i] = -1;
3970	independentOptions_[0] = 0;
3971	// But switch off slacks
3972	independentOptions_[1] = 512;
3973	// Substitute up to 3
3974	independentOptions_[2] = 3;
3975	}
3976	// Constructor when you really know what you are doing
3977	ClpSolve::ClpSolve(SolveType method, PresolveType presolveType,
3978	int numberPasses, int options[6],
3979	int extraInfo[6], int independentOptions[3])
3980	{
3981	method_ = method;
3982	presolveType_ = presolveType;
3983	numberPasses_ = numberPasses;
3984	int i;
3985	for (i = 0; i < 6; i++)
3986	options_[i] = options[i];
3987	options_[6] = 0;
3988	for (i = 0; i < 6; i++)
3989	extraInfo_[i] = extraInfo[i];
3990	extraInfo_[6] = 0;
3991	for (i = 0; i < 3; i++)
3992	independentOptions_[i] = independentOptions[i];
3993	}
3994
3995	// Copy constructor.
3996	ClpSolve::ClpSolve(const ClpSolve &rhs)
3997	{
3998	method_ = rhs.method_;
3999	presolveType_ = rhs.presolveType_;
4000	numberPasses_ = rhs.numberPasses_;
4001	int i;
4002	for (i = 0; i < 7; i++)
4003	options_[i] = rhs.options_[i];
4004	for (i = 0; i < 7; i++)
4005	extraInfo_[i] = rhs.extraInfo_[i];
4006	for (i = 0; i < 3; i++)
4007	independentOptions_[i] = rhs.independentOptions_[i];
4008	}
4009	// Assignment operator. This copies the data
4010	ClpSolve &
4011	ClpSolve::operator=(const ClpSolve &rhs)
4012	{
4013	if (this != &rhs) {
4014	method_ = rhs.method_;
4015	presolveType_ = rhs.presolveType_;
4016	numberPasses_ = rhs.numberPasses_;
4017	int i;
4018	for (i = 0; i < 7; i++)
4019	options_[i] = rhs.options_[i];
4020	for (i = 0; i < 7; i++)
4021	extraInfo_[i] = rhs.extraInfo_[i];
4022	for (i = 0; i < 3; i++)
4023	independentOptions_[i] = rhs.independentOptions_[i];
4024	}
4025	return *this;
4026	}
4027	// Destructor
4028	ClpSolve::~ClpSolve()
4029	{
4030	}
4031	// See header file for details
4032	void ClpSolve::setSpecialOption(int which, int value, int extraInfo)
4033	{
4034	options_[which] = value;
4035	extraInfo_[which] = extraInfo;
4036	}
4037	int ClpSolve::getSpecialOption(int which) const
4038	{
4039	return options_[which];
4040	}
4041
4042	// Solve types
4043	void ClpSolve::setSolveType(SolveType method, int /extraInfo/)
4044	{
4045	method_ = method;
4046	}
4047
4048	ClpSolve::SolveType
4049	ClpSolve::getSolveType()
4050	{
4051	return method_;
4052	}
4053
4054	// Presolve types
4055	void ClpSolve::setPresolveType(PresolveType amount, int extraInfo)
4056	{
4057	presolveType_ = amount;
4058	numberPasses_ = extraInfo;
4059	}
4060	ClpSolve::PresolveType
4061	ClpSolve::getPresolveType()
4062	{
4063	return presolveType_;
4064	}
4065	// Extra info for idiot (or sprint)
4066	int ClpSolve::getExtraInfo(int which) const
4067	{
4068	return extraInfo_[which];
4069	}
4070	int ClpSolve::getPresolvePasses() const
4071	{
4072	return numberPasses_;
4073	}
4074	/* Say to return at once if infeasible,
4075	default is to solve */
4076	void ClpSolve::setInfeasibleReturn(bool trueFalse)
4077	{
4078	independentOptions_[0] = trueFalse ? 1 : 0;
4079	}
4080	#include <string>
4081	// Generates code for above constructor
4082	void ClpSolve::generateCpp(FILE *fp)
4083	{
4084	std::string solveType[] = {
4085	"ClpSolve::useDual",
4086	"ClpSolve::usePrimal",
4087	"ClpSolve::usePrimalorSprint",
4088	"ClpSolve::useBarrier",
4089	"ClpSolve::useBarrierNoCross",
4090	"ClpSolve::automatic",
4091	"ClpSolve::notImplemented"
4092	};
4093	std::string presolveType[] = {
4094	"ClpSolve::presolveOn",
4095	"ClpSolve::presolveOff",
4096	"ClpSolve::presolveNumber",
4097	"ClpSolve::presolveNumberCost"
4098	};
4099	fprintf(fp, "3 ClpSolve::SolveType method = %s;\n", solveType[method_].c_str());
4100	fprintf(fp, "3 ClpSolve::PresolveType presolveType = %s;\n",
4101	presolveType[presolveType_].c_str());
4102	fprintf(fp, "3 int numberPasses = %d;\n", numberPasses_);
4103	fprintf(fp, "3 int options[] = {%d,%d,%d,%d,%d,%d};\n",
4104	options_[0], options_[1], options_[2],
4105	options_[3], options_[4], options_[5]);
4106	fprintf(fp, "3 int extraInfo[] = {%d,%d,%d,%d,%d,%d};\n",
4107	extraInfo_[0], extraInfo_[1], extraInfo_[2],
4108	extraInfo_[3], extraInfo_[4], extraInfo_[5]);
4109	fprintf(fp, "3 int independentOptions[] = {%d,%d,%d};\n",
4110	independentOptions_[0], independentOptions_[1], independentOptions_[2]);
4111	fprintf(fp, "3 ClpSolve clpSolve(method,presolveType,numberPasses,\n");
4112	fprintf(fp, "3 options,extraInfo,independentOptions);\n");
4113	}
4114	//#############################################################################
4115	#include "ClpNonLinearCost.hpp"
4116
4117	ClpSimplexProgress::ClpSimplexProgress()
4118	{
4119	int i;
4120	for (i = 0; i < CLP_PROGRESS; i++) {
4121	objective_[i] = COIN_DBL_MAX * 1.0e-50;
4122	infeasibility_[i] = -1.0; // set to an impossible value
4123	realInfeasibility_[i] = COIN_DBL_MAX * 1.0e-50;
4124	numberInfeasibilities_[i] = -1;
4125	iterationNumber_[i] = -1;
4126	}
4127	#ifdef CLP_PROGRESS_WEIGHT
4128	for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4129	objectiveWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4130	infeasibilityWeight_[i] = -1.0; // set to an impossible value
4131	realInfeasibilityWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4132	numberInfeasibilitiesWeight_[i] = -1;
4133	iterationNumberWeight_[i] = -1;
4134	}
4135	drop_ = 0.0;
4136	best_ = 0.0;
4137	#endif
4138	initialWeight_ = 0.0;
4139	for (i = 0; i < CLP_CYCLE; i++) {
4140	//obj_[i]=COIN_DBL_MAX*1.0e-50;
4141	in_[i] = -1;
4142	out_[i] = -1;
4143	way_[i] = 0;
4144	}
4145	numberTimes_ = 0;
4146	numberBadTimes_ = 0;
4147	numberReallyBadTimes_ = 0;
4148	numberTimesFlagged_ = 0;
4149	model_ = NULL;
4150	oddState_ = 0;
4151	}
4152
4153	//-----------------------------------------------------------------------------
4154
4155	ClpSimplexProgress::~ClpSimplexProgress()
4156	{
4157	}
4158	// Copy constructor.
4159	ClpSimplexProgress::ClpSimplexProgress(const ClpSimplexProgress &rhs)
4160	{
4161	int i;
4162	for (i = 0; i < CLP_PROGRESS; i++) {
4163	objective_[i] = rhs.objective_[i];
4164	infeasibility_[i] = rhs.infeasibility_[i];
4165	realInfeasibility_[i] = rhs.realInfeasibility_[i];
4166	numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
4167	iterationNumber_[i] = rhs.iterationNumber_[i];
4168	}
4169	#ifdef CLP_PROGRESS_WEIGHT
4170	for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4171	objectiveWeight_[i] = rhs.objectiveWeight_[i];
4172	infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
4173	realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
4174	numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
4175	iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
4176	}
4177	drop_ = rhs.drop_;
4178	best_ = rhs.best_;
4179	#endif
4180	initialWeight_ = rhs.initialWeight_;
4181	for (i = 0; i < CLP_CYCLE; i++) {
4182	//obj_[i]=rhs.obj_[i];
4183	in_[i] = rhs.in_[i];
4184	out_[i] = rhs.out_[i];
4185	way_[i] = rhs.way_[i];
4186	}
4187	numberTimes_ = rhs.numberTimes_;
4188	numberBadTimes_ = rhs.numberBadTimes_;
4189	numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
4190	numberTimesFlagged_ = rhs.numberTimesFlagged_;
4191	model_ = rhs.model_;
4192	oddState_ = rhs.oddState_;
4193	}
4194	// Copy constructor.from model
4195	ClpSimplexProgress::ClpSimplexProgress(ClpSimplex *model)
4196	{
4197	model_ = model;
4198	reset();
4199	initialWeight_ = 0.0;
4200	}
4201	// Fill from model
4202	void ClpSimplexProgress::fillFromModel(ClpSimplex *model)
4203	{
4204	model_ = model;
4205	reset();
4206	initialWeight_ = 0.0;
4207	}
4208	// Assignment operator. This copies the data
4209	ClpSimplexProgress &
4210	ClpSimplexProgress::operator=(const ClpSimplexProgress &rhs)
4211	{
4212	if (this != &rhs) {
4213	int i;
4214	for (i = 0; i < CLP_PROGRESS; i++) {
4215	objective_[i] = rhs.objective_[i];
4216	infeasibility_[i] = rhs.infeasibility_[i];
4217	realInfeasibility_[i] = rhs.realInfeasibility_[i];
4218	numberInfeasibilities_[i] = rhs.numberInfeasibilities_[i];
4219	iterationNumber_[i] = rhs.iterationNumber_[i];
4220	}
4221	#ifdef CLP_PROGRESS_WEIGHT
4222	for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4223	objectiveWeight_[i] = rhs.objectiveWeight_[i];
4224	infeasibilityWeight_[i] = rhs.infeasibilityWeight_[i];
4225	realInfeasibilityWeight_[i] = rhs.realInfeasibilityWeight_[i];
4226	numberInfeasibilitiesWeight_[i] = rhs.numberInfeasibilitiesWeight_[i];
4227	iterationNumberWeight_[i] = rhs.iterationNumberWeight_[i];
4228	}
4229	drop_ = rhs.drop_;
4230	best_ = rhs.best_;
4231	#endif
4232	initialWeight_ = rhs.initialWeight_;
4233	for (i = 0; i < CLP_CYCLE; i++) {
4234	//obj_[i]=rhs.obj_[i];
4235	in_[i] = rhs.in_[i];
4236	out_[i] = rhs.out_[i];
4237	way_[i] = rhs.way_[i];
4238	}
4239	numberTimes_ = rhs.numberTimes_;
4240	numberBadTimes_ = rhs.numberBadTimes_;
4241	numberReallyBadTimes_ = rhs.numberReallyBadTimes_;
4242	numberTimesFlagged_ = rhs.numberTimesFlagged_;
4243	model_ = rhs.model_;
4244	oddState_ = rhs.oddState_;
4245	}
4246	return *this;
4247	}
4248	// Seems to be something odd about exact comparison of doubles on linux
4249	static bool equalDouble(double value1, double value2)
4250	{
4251
4252	union {
4253	double d;
4254	int i[2];
4255	} v1, v2;
4256	v1.d = value1;
4257	v2.d = value2;
4258	if (sizeof(int) * 2 == sizeof(double))
4259	return (v1.i[0] == v2.i[0] && v1.i[1] == v2.i[1]);
4260	else
4261	return (v1.i[0] == v2.i[0]);
4262	}
4263	int ClpSimplexProgress::looping()
4264	{
4265	if (!model_)
4266	return -1;
4267	double objective;
4268	if (model_->algorithm() < 0) {
4269	objective = model_->rawObjectiveValue();
4270	objective -= model_->bestPossibleImprovement();
4271	} else {
4272	objective = model_->nonLinearCost()->feasibleReportCost();
4273	}
4274	double infeasibility;
4275	double realInfeasibility = 0.0;
4276	int numberInfeasibilities;
4277	int iterationNumber = model_->numberIterations();
4278	//numberTimesFlagged_ = 0;
4279	if (model_->algorithm() < 0) {
4280	// dual
4281	infeasibility = model_->sumPrimalInfeasibilities();
4282	numberInfeasibilities = model_->numberPrimalInfeasibilities();
4283	} else {
4284	//primal
4285	infeasibility = model_->sumDualInfeasibilities();
4286	realInfeasibility = model_->nonLinearCost()->sumInfeasibilities();
4287	numberInfeasibilities = model_->numberDualInfeasibilities();
4288	}
4289	int i;
4290	int numberMatched = 0;
4291	int matched = 0;
4292	int nsame = 0;
4293	for (i = 0; i < CLP_PROGRESS; i++) {
4294	bool matchedOnObjective = equalDouble(objective, objective_[i]);
4295	bool matchedOnInfeasibility = equalDouble(infeasibility, infeasibility_[i]);
4296	bool matchedOnInfeasibilities = (numberInfeasibilities == numberInfeasibilities_[i]);
4297
4298	if (matchedOnObjective && matchedOnInfeasibility && matchedOnInfeasibilities) {
4299	matched \|= (1 << i);
4300	// Check not same iteration
4301	if (iterationNumber != iterationNumber_[i]) {
4302	numberMatched++;
4303	// here mainly to get over compiler bug?
4304	if (model_->messageHandler()->logLevel() > 10)
4305	printf("%d %d %d %d %d loop check\n", i, numberMatched,
4306	matchedOnObjective, matchedOnInfeasibility,
4307	matchedOnInfeasibilities);
4308	} else {
4309	// stuck but code should notice
4310	nsame++;
4311	}
4312	}
4313	if (i) {
4314	objective_[i - 1] = objective_[i];
4315	infeasibility_[i - 1] = infeasibility_[i];
4316	realInfeasibility_[i - 1] = realInfeasibility_[i];
4317	numberInfeasibilities_[i - 1] = numberInfeasibilities_[i];
4318	iterationNumber_[i - 1] = iterationNumber_[i];
4319	}
4320	}
4321	objective_[CLP_PROGRESS - 1] = objective;
4322	infeasibility_[CLP_PROGRESS - 1] = infeasibility;
4323	realInfeasibility_[CLP_PROGRESS - 1] = realInfeasibility;
4324	numberInfeasibilities_[CLP_PROGRESS - 1] = numberInfeasibilities;
4325	iterationNumber_[CLP_PROGRESS - 1] = iterationNumber;
4326	if (nsame == CLP_PROGRESS)
4327	numberMatched = CLP_PROGRESS; // really stuck
4328	if (model_->progressFlag())
4329	numberMatched = 0;
4330	numberTimes_++;
4331	if (numberTimes_ < 10)
4332	numberMatched = 0;
4333	// skip if just last time as may be checking something
4334	if (matched == (1 << (CLP_PROGRESS - 1)))
4335	numberMatched = 0;
4336	if (numberMatched && model_->clpMatrix()->type() < 15) {
4337	model_->messageHandler()->message(CLP_POSSIBLELOOP, model_->messages())
4338	<< numberMatched
4339	<< matched
4340	<< numberTimes_
4341	<< CoinMessageEol;
4342	numberBadTimes_++;
4343	if (numberBadTimes_ < 10) {
4344	// make factorize every iteration
4345	model_->forceFactorization(1);
4346	if (numberBadTimes_ < 2) {
4347	startCheck(); // clear other loop check
4348	if (model_->algorithm() < 0) {
4349	// dual - change tolerance
4350	model_->setCurrentDualTolerance(model_->currentDualTolerance() * 1.05);
4351	// if infeasible increase dual bound
4352	if (model_->dualBound() < 1.0e17) {
4353	model_->setDualBound(model_->dualBound() * 1.1);
4354	static_cast< ClpSimplexDual * >(model_)->resetFakeBounds(0);
4355	}
4356	} else {
4357	// primal - change tolerance
4358	if (numberBadTimes_ > 3)
4359	model_->setCurrentPrimalTolerance(model_->currentPrimalTolerance() * 1.05);
4360	// if infeasible increase infeasibility cost
4361	if (model_->nonLinearCost()->numberInfeasibilities() && model_->infeasibilityCost() < 1.0e17) {
4362	model_->setInfeasibilityCost(model_->infeasibilityCost() * 1.1);
4363	}
4364	}
4365	} else {
4366	// flag
4367	int iSequence;
4368	if (model_->algorithm() < 0) {
4369	// dual
4370	if (model_->dualBound() > 1.0e14)
4371	model_->setDualBound(1.0e14);
4372	iSequence = in_[CLP_CYCLE - 1];
4373	} else {
4374	// primal
4375	//if (model_->infeasibilityCost() > 1.0e14)
4376	// model_->setInfeasibilityCost(1.0e14);
4377	iSequence = out_[CLP_CYCLE - 1];
4378	}
4379	if (iSequence >= 0) {
4380	char x = model_->isColumn(iSequence) ? 'C' : 'R';
4381	if (model_->messageHandler()->logLevel() >= 63)
4382	model_->messageHandler()->message(CLP_SIMPLEX_FLAG, model_->messages())
4383	<< x << model_->sequenceWithin(iSequence)
4384	<< CoinMessageEol;
4385	// if Gub then needs to be sequenceIn_
4386	int save = model_->sequenceIn();
4387	model_->setSequenceIn(iSequence);
4388	model_->setFlagged(iSequence);
4389	model_->setSequenceIn(save);
4390	//printf("flagging %d from loop\n",iSequence);
4391	startCheck();
4392	} else {
4393	// Give up
4394	if (model_->messageHandler()->logLevel() >= 63)
4395	printf("***** All flagged?\n");
4396	return 4;
4397	}
4398	// reset
4399	numberBadTimes_ = 2;
4400	}
4401	return -2;
4402	} else {
4403	// look at solution and maybe declare victory
4404	if (infeasibility < 1.0e-4) {
4405	return 0;
4406	} else {
4407	model_->messageHandler()->message(CLP_LOOP, model_->messages())
4408	<< CoinMessageEol;
4409	#ifndef NDEBUG
4410	printf("debug loop ClpSimplex A\n");
4411	abort();
4412	#endif
4413	return 3;
4414	}
4415	}
4416	}
4417	return -1;
4418	}
4419	// Resets as much as possible
4420	void ClpSimplexProgress::reset()
4421	{
4422	int i;
4423	for (i = 0; i < CLP_PROGRESS; i++) {
4424	if (model_->algorithm() >= 0)
4425	objective_[i] = COIN_DBL_MAX * 1.0e-50;
4426	else
4427	objective_[i] = -COIN_DBL_MAX * 1.0e-50;
4428	infeasibility_[i] = -1.0; // set to an impossible value
4429	realInfeasibility_[i] = COIN_DBL_MAX * 1.0e-50;
4430	numberInfeasibilities_[i] = -1;
4431	iterationNumber_[i] = -1;
4432	}
4433	#ifdef CLP_PROGRESS_WEIGHT
4434	for (i = 0; i < CLP_PROGRESS_WEIGHT; i++) {
4435	objectiveWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4436	infeasibilityWeight_[i] = -1.0; // set to an impossible value
4437	realInfeasibilityWeight_[i] = COIN_DBL_MAX * 1.0e-50;
4438	numberInfeasibilitiesWeight_[i] = -1;
4439	iterationNumberWeight_[i] = -1;
4440	}
4441	drop_ = 0.0;
4442	best_ = 0.0;
4443	#endif
4444	for (i = 0; i < CLP_CYCLE; i++) {
4445	//obj_[i]=COIN_DBL_MAX*1.0e-50;
4446	in_[i] = -1;
4447	out_[i] = -1;
4448	way_[i] = 0;
4449	}
4450	numberTimes_ = 0;
4451	numberBadTimes_ = 0;
4452	numberReallyBadTimes_ = 0;
4453	numberTimesFlagged_ = 0;
4454	oddState_ = 0;
4455	}
4456	// Returns previous objective (if -1) - current if (0)
4457	double
4458	ClpSimplexProgress::lastObjective(int back) const
4459	{
4460	return objective_[CLP_PROGRESS - 1 - back];
4461	}
4462	// Returns previous infeasibility (if -1) - current if (0)
4463	double
4464	ClpSimplexProgress::lastInfeasibility(int back) const
4465	{
4466	return realInfeasibility_[CLP_PROGRESS - 1 - back];
4467	}
4468	// Sets real primal infeasibility
4469	void ClpSimplexProgress::setInfeasibility(double value)
4470	{
4471	for (int i = 1; i < CLP_PROGRESS; i++)
4472	realInfeasibility_[i - 1] = realInfeasibility_[i];
4473	realInfeasibility_[CLP_PROGRESS - 1] = value;
4474	}
4475	// Returns number of primal infeasibilities (if -1) - current if (0)
4476	int ClpSimplexProgress::numberInfeasibilities(int back) const
4477	{
4478	return numberInfeasibilities_[CLP_PROGRESS - 1 - back];
4479	}
4480	// Modify objective e.g. if dual infeasible in dual
4481	void ClpSimplexProgress::modifyObjective(double value)
4482	{
4483	objective_[CLP_PROGRESS - 1] = value;
4484	}
4485	// Returns previous iteration number (if -1) - current if (0)
4486	int ClpSimplexProgress::lastIterationNumber(int back) const
4487	{
4488	return iterationNumber_[CLP_PROGRESS - 1 - back];
4489	}
4490	// clears iteration numbers (to switch off panic)
4491	void ClpSimplexProgress::clearIterationNumbers()
4492	{
4493	for (int i = 0; i < CLP_PROGRESS; i++)
4494	iterationNumber_[i] = -1;
4495	}
4496	// Start check at beginning of whileIterating
4497	void ClpSimplexProgress::startCheck()
4498	{
4499	int i;
4500	for (i = 0; i < CLP_CYCLE; i++) {
4501	//obj_[i]=COIN_DBL_MAX*1.0e-50;
4502	in_[i] = -1;
4503	out_[i] = -1;
4504	way_[i] = 0;
4505	}
4506	}
4507	// Returns cycle length in whileIterating
4508	int ClpSimplexProgress::cycle(int in, int out, int wayIn, int wayOut)
4509	{
4510	int i;
4511	#if 0
4512	if (model_->numberIterations() > 206571) {
4513	printf("in %d out %d\n", in, out);
4514	for (i = 0; i < CLP_CYCLE; i++)
4515	printf("cy %d in %d out %d\n", i, in_[i], out_[i]);
4516	}
4517	#endif
4518	int matched = 0;
4519	// first see if in matches any out
4520	for (i = 1; i < CLP_CYCLE; i++) {
4521	if (in == out_[i]) {
4522	// even if flip then suspicious
4523	matched = -1;
4524	break;
4525	}
4526	}
4527	#if 0
4528	if (!matched \|\| in_[0] < 0) {
4529	// can't be cycle
4530	for (i = 0; i < CLP_CYCLE - 1; i++) {
4531	//obj_[i]=obj_[i+1];
4532	in_[i] = in_[i+1];
4533	out_[i] = out_[i+1];
4534	way_[i] = way_[i+1];
4535	}
4536	} else {
4537	// possible cycle
4538	matched = 0;
4539	for (i = 0; i < CLP_CYCLE - 1; i++) {
4540	int k;
4541	char wayThis = way_[i];
4542	int inThis = in_[i];
4543	int outThis = out_[i];
4544	//double objThis = obj_[i];
4545	for(k = i + 1; k < CLP_CYCLE; k++) {
4546	if (inThis == in_[k] && outThis == out_[k] && wayThis == way_[k]) {
4547	int distance = k - i;
4548	if (k + distance < CLP_CYCLE) {
4549	// See if repeats
4550	int j = k + distance;
4551	if (inThis == in_[j] && outThis == out_[j] && wayThis == way_[j]) {
4552	matched = distance;
4553	break;
4554	}
4555	} else {
4556	matched = distance;
4557	break;
4558	}
4559	}
4560	}
4561	//obj_[i]=obj_[i+1];
4562	in_[i] = in_[i+1];
4563	out_[i] = out_[i+1];
4564	way_[i] = way_[i+1];
4565	}
4566	}
4567	#else
4568	if (matched && in_[0] >= 0) {
4569	// possible cycle - only check [0] against all
4570	matched = 0;
4571	int nMatched = 0;
4572	char way0 = way_[0];
4573	int in0 = in_[0];
4574	int out0 = out_[0];
4575	//double obj0 = obj_[i];
4576	for (int k = 1; k < CLP_CYCLE - 4; k++) {
4577	if (in0 == in_[k] && out0 == out_[k] && way0 == way_[k]) {
4578	nMatched++;
4579	// See if repeats
4580	int end = CLP_CYCLE - k;
4581	int j;
4582	for (j = 1; j < end; j++) {
4583	if (in_[j + k] != in_[j] \|\| out_[j + k] != out_[j] \|\| way_[j + k] != way_[j])
4584	break;
4585	}
4586	if (j == end) {
4587	matched = k;
4588	break;
4589	}
4590	}
4591	}
4592	// If three times then that is too much even if not regular
4593	if (matched <= 0 && nMatched > 1)
4594	matched = 100;
4595	}
4596	for (i = 0; i < CLP_CYCLE - 1; i++) {
4597	//obj_[i]=obj_[i+1];
4598	in_[i] = in_[i + 1];
4599	out_[i] = out_[i + 1];
4600	way_[i] = way_[i + 1];
4601	}
4602	#endif
4603	int way = 1 - wayIn + 4 * (1 - wayOut);
4604	//obj_[i]=model_->objectiveValue();
4605	in_[CLP_CYCLE - 1] = in;
4606	out_[CLP_CYCLE - 1] = out;
4607	way_[CLP_CYCLE - 1] = static_cast< char >(way);
4608	return matched;
4609	}
4610	#include "CoinStructuredModel.hpp"
4611	// Solve using structure of model and maybe in parallel
4612	int ClpSimplex::solve(CoinStructuredModel *model)
4613	{
4614	// analyze structure
4615	int numberRowBlocks = model->numberRowBlocks();
4616	int numberColumnBlocks = model->numberColumnBlocks();
4617	int numberElementBlocks = model->numberElementBlocks();
4618	if (numberElementBlocks == 1) {
4619	loadProblem(*model, false);
4620	return dual();
4621	}
4622	// For now just get top level structure
4623	CoinModelBlockInfo *blockInfo = new CoinModelBlockInfo[numberElementBlocks];
4624	for (int i = 0; i < numberElementBlocks; i++) {
4625	CoinStructuredModel subModel = dynamic_cast< CoinStructuredModel >(model->block(i));
4626	CoinModel *thisBlock;
4627	if (subModel) {
4628	thisBlock = subModel->coinModelBlock(blockInfo[i]);
4629	model->setCoinModel(thisBlock, i);
4630	} else {
4631	thisBlock = dynamic_cast< CoinModel * >(model->block(i));
4632	assert(thisBlock);
4633	// just fill in info
4634	CoinModelBlockInfo info = CoinModelBlockInfo();
4635	int whatsSet = thisBlock->whatIsSet();
4636	info.matrix = static_cast< char >(((whatsSet & 1) != 0) ? 1 : 0);
4637	info.rhs = static_cast< char >(((whatsSet & 2) != 0) ? 1 : 0);
4638	info.rowName = static_cast< char >(((whatsSet & 4) != 0) ? 1 : 0);
4639	info.integer = static_cast< char >(((whatsSet & 32) != 0) ? 1 : 0);
4640	info.bounds = static_cast< char >(((whatsSet & 8) != 0) ? 1 : 0);
4641	info.columnName = static_cast< char >(((whatsSet & 16) != 0) ? 1 : 0);
4642	// Which block
4643	int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
4644	info.rowBlock = iRowBlock;
4645	int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
4646	info.columnBlock = iColumnBlock;
4647	blockInfo[i] = info;
4648	}
4649	}
4650	int *rowCounts = new int[numberRowBlocks];
4651	CoinZeroN(rowCounts, numberRowBlocks);
4652	int *columnCounts = new int[numberColumnBlocks + 1];
4653	CoinZeroN(columnCounts, numberColumnBlocks);
4654	int decomposeType = 0;
4655	for (int i = 0; i < numberElementBlocks; i++) {
4656	int iRowBlock = blockInfo[i].rowBlock;
4657	int iColumnBlock = blockInfo[i].columnBlock;
4658	rowCounts[iRowBlock]++;
4659	columnCounts[iColumnBlock]++;
4660	}
4661	if (numberRowBlocks == numberColumnBlocks \|\| numberRowBlocks == numberColumnBlocks + 1) {
4662	// could be Dantzig-Wolfe
4663	int numberG1 = 0;
4664	for (int i = 0; i < numberRowBlocks; i++) {
4665	if (rowCounts[i] > 1)
4666	numberG1++;
4667	}
4668	bool masterColumns = (numberColumnBlocks == numberRowBlocks);
4669	if ((masterColumns && numberElementBlocks == 2 * numberRowBlocks - 1)
4670	\|\| (!masterColumns && numberElementBlocks == 2 * numberRowBlocks)) {
4671	if (numberG1 < 2)
4672	decomposeType = 1;
4673	}
4674	}
4675	if (!decomposeType && (numberRowBlocks == numberColumnBlocks \|\| numberRowBlocks == numberColumnBlocks - 1)) {
4676	// could be Benders
4677	int numberG1 = 0;
4678	for (int i = 0; i < numberColumnBlocks; i++) {
4679	if (columnCounts[i] > 1)
4680	numberG1++;
4681	}
4682	bool masterRows = (numberColumnBlocks == numberRowBlocks);
4683	if ((masterRows && numberElementBlocks == 2 * numberColumnBlocks - 1)
4684	\|\| (!masterRows && numberElementBlocks == 2 * numberColumnBlocks)) {
4685	if (numberG1 < 2)
4686	decomposeType = 2;
4687	}
4688	}
4689	delete[] rowCounts;
4690	delete[] columnCounts;
4691	delete[] blockInfo;
4692	// decide what to do
4693	ClpSolve options;
4694	options.setIndependentOption(2, 100);
4695	switch (decomposeType) {
4696	// No good
4697	case 0:
4698	loadProblem(*model, false);
4699	return dual();
4700	// DW
4701	case 1:
4702	return solveDW(model, options);
4703	// Benders
4704	case 2:
4705	return solveBenders(model, options);
4706	}
4707	return 0; // to stop compiler warning
4708	}
4709	/* This loads a model from a CoinStructuredModel object - returns number of errors.
4710	If originalOrder then keep to order stored in blocks,
4711	otherwise first column/rows correspond to first block - etc.
4712	If keepSolution true and size is same as current then
4713	keeps current status and solution
4714	*/
4715	int ClpSimplex::loadProblem(CoinStructuredModel &coinModel,
4716	bool originalOrder,
4717	bool keepSolution)
4718	{
4719	unsigned char *status = NULL;
4720	double *psol = NULL;
4721	double *dsol = NULL;
4722	int numberRows = coinModel.numberRows();
4723	int numberColumns = coinModel.numberColumns();
4724	int numberRowBlocks = coinModel.numberRowBlocks();
4725	int numberColumnBlocks = coinModel.numberColumnBlocks();
4726	int numberElementBlocks = coinModel.numberElementBlocks();
4727	if (status_ && numberRows_ && numberRows_ == numberRows && numberColumns_ == numberColumns && keepSolution) {
4728	status = new unsigned char[numberRows_ + numberColumns_];
4729	CoinMemcpyN(status_, numberRows_ + numberColumns_, status);
4730	psol = new double[numberRows_ + numberColumns_];
4731	CoinMemcpyN(columnActivity_, numberColumns_, psol);
4732	CoinMemcpyN(rowActivity_, numberRows_, psol + numberColumns_);
4733	dsol = new double[numberRows_ + numberColumns_];
4734	CoinMemcpyN(reducedCost_, numberColumns_, dsol);
4735	CoinMemcpyN(dual_, numberRows_, dsol + numberColumns_);
4736	}
4737	int returnCode = 0;
4738	double *rowLower = new double[numberRows];
4739	double *rowUpper = new double[numberRows];
4740	double *columnLower = new double[numberColumns];
4741	double *columnUpper = new double[numberColumns];
4742	double *objective = new double[numberColumns];
4743	int *integerType = new int[numberColumns];
4744	CoinBigIndex numberElements = 0;
4745	// Bases for blocks
4746	int *rowBase = new int[numberRowBlocks];
4747	CoinFillN(rowBase, numberRowBlocks, -1);
4748	// And row to put it
4749	int *whichRow = new int[numberRows + numberRowBlocks];
4750	int *columnBase = new int[numberColumnBlocks];
4751	CoinFillN(columnBase, numberColumnBlocks, -1);
4752	// And column to put it
4753	int *whichColumn = new int[numberColumns + numberColumnBlocks];
4754	for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4755	CoinModel *block = coinModel.coinBlock(iBlock);
4756	numberElements += block->numberElements();
4757	//and set up elements etc
4758	double *associated = block->associatedArray();
4759	// If strings then do copies
4760	if (block->stringsExist())
4761	returnCode += block->createArrays(rowLower, rowUpper, columnLower, columnUpper,
4762	objective, integerType, associated);
4763	const CoinModelBlockInfo &info = coinModel.blockType(iBlock);
4764	int iRowBlock = info.rowBlock;
4765	int iColumnBlock = info.columnBlock;
4766	if (rowBase[iRowBlock] < 0) {
4767	rowBase[iRowBlock] = block->numberRows();
4768	// Save block number
4769	whichRow[numberRows + iRowBlock] = iBlock;
4770	} else {
4771	assert(rowBase[iRowBlock] == block->numberRows());
4772	}
4773	if (columnBase[iColumnBlock] < 0) {
4774	columnBase[iColumnBlock] = block->numberColumns();
4775	// Save block number
4776	whichColumn[numberColumns + iColumnBlock] = iBlock;
4777	} else {
4778	assert(columnBase[iColumnBlock] == block->numberColumns());
4779	}
4780	}
4781	// Fill arrays with defaults
4782	CoinFillN(rowLower, numberRows, -COIN_DBL_MAX);
4783	CoinFillN(rowUpper, numberRows, COIN_DBL_MAX);
4784	CoinFillN(columnLower, numberColumns, 0.0);
4785	CoinFillN(columnUpper, numberColumns, COIN_DBL_MAX);
4786	CoinFillN(objective, numberColumns, 0.0);
4787	CoinFillN(integerType, numberColumns, 0);
4788	int n = 0;
4789	for (int iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
4790	int k = rowBase[iBlock];
4791	rowBase[iBlock] = n;
4792	assert(k >= 0);
4793	// block number
4794	int jBlock = whichRow[numberRows + iBlock];
4795	if (originalOrder) {
4796	memcpy(whichRow + n, coinModel.coinBlock(jBlock)->originalRows(), k * sizeof(int));
4797	} else {
4798	CoinIotaN(whichRow + n, k, n);
4799	}
4800	n += k;
4801	}
4802	assert(n == numberRows);
4803	n = 0;
4804	for (int iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
4805	int k = columnBase[iBlock];
4806	columnBase[iBlock] = n;
4807	assert(k >= 0);
4808	if (k) {
4809	// block number
4810	int jBlock = whichColumn[numberColumns + iBlock];
4811	if (originalOrder) {
4812	memcpy(whichColumn + n, coinModel.coinBlock(jBlock)->originalColumns(),
4813	k * sizeof(int));
4814	} else {
4815	CoinIotaN(whichColumn + n, k, n);
4816	}
4817	n += k;
4818	}
4819	}
4820	assert(n == numberColumns);
4821	bool gotIntegers = false;
4822	for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4823	CoinModel *block = coinModel.coinBlock(iBlock);
4824	const CoinModelBlockInfo &info = coinModel.blockType(iBlock);
4825	int iRowBlock = info.rowBlock;
4826	int iRowBase = rowBase[iRowBlock];
4827	int iColumnBlock = info.columnBlock;
4828	int iColumnBase = columnBase[iColumnBlock];
4829	if (info.rhs) {
4830	int nRows = block->numberRows();
4831	const double *lower = block->rowLowerArray();
4832	const double *upper = block->rowUpperArray();
4833	for (int i = 0; i < nRows; i++) {
4834	int put = whichRow[i + iRowBase];
4835	rowLower[put] = lower[i];
4836	rowUpper[put] = upper[i];
4837	}
4838	}
4839	if (info.bounds) {
4840	int nColumns = block->numberColumns();
4841	const double *lower = block->columnLowerArray();
4842	const double *upper = block->columnUpperArray();
4843	const double *obj = block->objectiveArray();
4844	for (int i = 0; i < nColumns; i++) {
4845	int put = whichColumn[i + iColumnBase];
4846	columnLower[put] = lower[i];
4847	columnUpper[put] = upper[i];
4848	objective[put] = obj[i];
4849	}
4850	}
4851	if (info.integer) {
4852	gotIntegers = true;
4853	int nColumns = block->numberColumns();
4854	const int *type = block->integerTypeArray();
4855	for (int i = 0; i < nColumns; i++) {
4856	int put = whichColumn[i + iColumnBase];
4857	integerType[put] = type[i];
4858	}
4859	}
4860	}
4861	gutsOfLoadModel(numberRows, numberColumns,
4862	columnLower, columnUpper, objective, rowLower, rowUpper, NULL);
4863	delete[] rowLower;
4864	delete[] rowUpper;
4865	delete[] columnLower;
4866	delete[] columnUpper;
4867	delete[] objective;
4868	// Do integers if wanted
4869	if (gotIntegers) {
4870	for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
4871	if (integerType[iColumn])
4872	setInteger(iColumn);
4873	}
4874	}
4875	delete[] integerType;
4876	setObjectiveOffset(coinModel.objectiveOffset());
4877	// Space for elements
4878	int *row = new int[numberElements];
4879	int *column = new int[numberElements];
4880	double *element = new double[numberElements];
4881	numberElements = 0;
4882	for (int iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
4883	CoinModel *block = coinModel.coinBlock(iBlock);
4884	const CoinModelBlockInfo &info = coinModel.blockType(iBlock);
4885	int iRowBlock = info.rowBlock;
4886	int iRowBase = rowBase[iRowBlock];
4887	int iColumnBlock = info.columnBlock;
4888	int iColumnBase = columnBase[iColumnBlock];
4889	if (info.rowName) {
4890	int numberItems = block->rowNames()->numberItems();
4891	assert(block->numberRows() >= numberItems);
4892	if (numberItems) {
4893	const char const rowNames = block->rowNames()->names();
4894	for (int i = 0; i < numberItems; i++) {
4895	int put = whichRow[i + iRowBase];
4896	std::string name = rowNames[i];
4897	setRowName(put, name);
4898	}
4899	}
4900	}
4901	if (info.columnName) {
4902	int numberItems = block->columnNames()->numberItems();
4903	assert(block->numberColumns() >= numberItems);
4904	if (numberItems) {
4905	const char const columnNames = block->columnNames()->names();
4906	for (int i = 0; i < numberItems; i++) {
4907	int put = whichColumn[i + iColumnBase];
4908	std::string name = columnNames[i];
4909	setColumnName(put, name);
4910	}
4911	}
4912	}
4913	if (info.matrix) {
4914	CoinPackedMatrix matrix2;
4915	const CoinPackedMatrix *matrix = block->packedMatrix();
4916	if (!matrix) {
4917	double *associated = block->associatedArray();
4918	block->createPackedMatrix(matrix2, associated);
4919	matrix = &matrix2;
4920	}
4921	// get matrix data pointers
4922	const int *row2 = matrix->getIndices();
4923	const CoinBigIndex *columnStart = matrix->getVectorStarts();
4924	const double *elementByColumn = matrix->getElements();
4925	const int *columnLength = matrix->getVectorLengths();
4926	int n = matrix->getNumCols();
4927	assert(matrix->isColOrdered());
4928	for (int iColumn = 0; iColumn < n; iColumn++) {
4929	CoinBigIndex j;
4930	int jColumn = whichColumn[iColumn + iColumnBase];
4931	for (j = columnStart[iColumn];
4932	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4933	row[numberElements] = whichRow[row2[j] + iRowBase];
4934	column[numberElements] = jColumn;
4935	element[numberElements++] = elementByColumn[j];
4936	}
4937	}
4938	}
4939	}
4940	delete[] whichRow;
4941	delete[] whichColumn;
4942	delete[] rowBase;
4943	delete[] columnBase;
4944	CoinPackedMatrix *matrix = new CoinPackedMatrix(true, row, column, element, numberElements);
4945	matrix_ = new ClpPackedMatrix(matrix);
4946	matrix_->setDimensions(numberRows, numberColumns);
4947	delete[] row;
4948	delete[] column;
4949	delete[] element;
4950	createStatus();
4951	if (status) {
4952	// copy back
4953	CoinMemcpyN(status, numberRows_ + numberColumns_, status_);
4954	CoinMemcpyN(psol, numberColumns_, columnActivity_);
4955	CoinMemcpyN(psol + numberColumns_, numberRows_, rowActivity_);
4956	CoinMemcpyN(dsol, numberColumns_, reducedCost_);
4957	CoinMemcpyN(dsol + numberColumns_, numberRows_, dual_);
4958	delete[] status;
4959	delete[] psol;
4960	delete[] dsol;
4961	}
4962	optimizationDirection_ = coinModel.optimizationDirection();
4963	return returnCode;
4964	}
4965	/* If input negative scales objective so maximum <= -value
4966	and returns scale factor used. If positive unscales and also
4967	redoes dual stuff
4968	*/
4969	double
4970	ClpSimplex::scaleObjective(double value)
4971	{
4972	double *obj = objective();
4973	double largest = 0.0;
4974	if (value < 0.0) {
4975	value = -value;
4976	for (int i = 0; i < numberColumns_; i++) {
4977	largest = CoinMax(largest, fabs(obj[i]));
4978	}
4979	if (largest > value) {
4980	double scaleFactor = value / largest;
4981	for (int i = 0; i < numberColumns_; i++) {
4982	obj[i] *= scaleFactor;
4983	reducedCost_[i] *= scaleFactor;
4984	}
4985	for (int i = 0; i < numberRows_; i++) {
4986	dual_[i] *= scaleFactor;
4987	}
4988	largest /= value;
4989	} else {
4990	// no need
4991	largest = 1.0;
4992	}
4993	} else {
4994	// at end
4995	if (value != 1.0) {
4996	for (int i = 0; i < numberColumns_; i++) {
4997	obj[i] *= value;
4998	reducedCost_[i] *= value;
4999	}
5000	for (int i = 0; i < numberRows_; i++) {
5001	dual_[i] *= value;
5002	}
5003	computeObjectiveValue();
5004	}
5005	}
5006	return largest;
5007	}
5008	#if defined(ABC_INHERIT) \|\| defined(CBC_THREAD) \|\| defined(THREADS_IN_ANALYZE)
5009	void clp_parallelManager(void stuff)
5010	{
5011	CoinPthreadStuff driver = reinterpret_cast< CoinPthreadStuff >(stuff);
5012	int whichThread = driver->whichThread();
5013	CoinThreadInfo *threadInfo = driver->threadInfoPointer(whichThread);
5014	threadInfo->status = -1;
5015	int *which = threadInfo->stuff;
5016	#ifdef PTHREAD_BARRIER_SERIAL_THREAD
5017	pthread_barrier_wait(driver->barrierPointer());
5018	#endif
5019	#if 0
5020	int status=-1;
5021	while (status!=100)
5022	status=timedWait(driver,1000,2);
5023	pthread_cond_signal(driver->conditionPointer(1));
5024	pthread_mutex_unlock(driver->mutexPointer(1,whichThread));
5025	#endif
5026	// so now mutex_ is locked
5027	int whichLocked = 0;
5028	while (true) {
5029	pthread_mutex_t *mutexPointer = driver->mutexPointer(whichLocked, whichThread);
5030	// wait
5031	//printf("Child waiting for %d - status %d %d %d\n",
5032	// whichLocked,lockedX[0],lockedX[1],lockedX[2]);
5033	#ifdef DETAIL_THREAD
5034	printf("thread %d about to lock mutex %d\n", whichThread, whichLocked);
5035	#endif
5036	pthread_mutex_lock(mutexPointer);
5037	whichLocked++;
5038	if (whichLocked == 3)
5039	whichLocked = 0;
5040	int unLock = whichLocked + 1;
5041	if (unLock == 3)
5042	unLock = 0;
5043	//printf("child pointer %x status %d\n",threadInfo,threadInfo->status);
5044	assert(threadInfo->status >= 0);
5045	if (threadInfo->status == 1000)
5046	pthread_exit(NULL);
5047	int type = threadInfo->status;
5048	int &returnCode = which[0];
5049	int iPass = which[1];
5050	ClpSimplex clpSimplex = reinterpret_cast< ClpSimplex >(threadInfo->extraInfo);
5051	//CoinIndexedVector * array;
5052	//double dummy;
5053	switch (type) {
5054	// dummy
5055	case 0:
5056	break;
5057	case 1:
5058	if (!clpSimplex->problemStatus() \|\| !iPass)
5059	returnCode = clpSimplex->dual();
5060	else
5061	returnCode = clpSimplex->primal();
5062	break;
5063	case 100:
5064	// initialization
5065	break;
5066	}
5067	threadInfo->status = -1;
5068	#ifdef DETAIL_THREAD
5069	printf("thread %d about to unlock mutex %d\n", whichThread, unLock);
5070	#endif
5071	pthread_mutex_unlock(driver->mutexPointer(unLock, whichThread));
5072	}
5073	}
5074	#endif
5075	// Solve using Dantzig-Wolfe decomposition and maybe in parallel
5076	int ClpSimplex::solveDW(CoinStructuredModel *model, ClpSolve &options)
5077	{
5078	double time1 = CoinCpuTime();
5079	int numberColumns = model->numberColumns();
5080	int numberRowBlocks = model->numberRowBlocks();
5081	int numberColumnBlocks = model->numberColumnBlocks();
5082	int numberElementBlocks = model->numberElementBlocks();
5083	// We already have top level structure
5084	CoinModelBlockInfo *blockInfo = new CoinModelBlockInfo[numberElementBlocks];
5085	for (int i = 0; i < numberElementBlocks; i++) {
5086	CoinModel *thisBlock = model->coinBlock(i);
5087	assert(thisBlock);
5088	// just fill in info
5089	CoinModelBlockInfo info = CoinModelBlockInfo();
5090	int whatsSet = thisBlock->whatIsSet();
5091	info.matrix = static_cast< char >(((whatsSet & 1) != 0) ? 1 : 0);
5092	info.rhs = static_cast< char >(((whatsSet & 2) != 0) ? 1 : 0);
5093	info.rowName = static_cast< char >(((whatsSet & 4) != 0) ? 1 : 0);
5094	info.integer = static_cast< char >(((whatsSet & 32) != 0) ? 1 : 0);
5095	info.bounds = static_cast< char >(((whatsSet & 8) != 0) ? 1 : 0);
5096	info.columnName = static_cast< char >(((whatsSet & 16) != 0) ? 1 : 0);
5097	// Which block
5098	int iRowBlock = model->rowBlock(thisBlock->getRowBlock());
5099	info.rowBlock = iRowBlock;
5100	int iColumnBlock = model->columnBlock(thisBlock->getColumnBlock());
5101	info.columnBlock = iColumnBlock;
5102	blockInfo[i] = info;
5103	}
5104	// make up problems
5105	int numberBlocks = numberRowBlocks - 1;
5106	// Find master rows and columns
5107	int *rowCounts = new int[numberRowBlocks];
5108	CoinZeroN(rowCounts, numberRowBlocks);
5109	int *columnCounts = new int[numberColumnBlocks + 1];
5110	CoinZeroN(columnCounts, numberColumnBlocks);
5111	int iBlock;
5112	for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5113	int iRowBlock = blockInfo[iBlock].rowBlock;
5114	rowCounts[iRowBlock]++;
5115	int iColumnBlock = blockInfo[iBlock].columnBlock;
5116	columnCounts[iColumnBlock]++;
5117	}
5118	int *whichBlock = new int[numberElementBlocks];
5119	int masterRowBlock = -1;
5120	for (iBlock = 0; iBlock < numberRowBlocks; iBlock++) {
5121	if (rowCounts[iBlock] > 1) {
5122	if (masterRowBlock == -1) {
5123	masterRowBlock = iBlock;
5124	} else {
5125	// Can't decode
5126	masterRowBlock = -2;
5127	break;
5128	}
5129	}
5130	}
5131	int masterColumnBlock = -1;
5132	int kBlock = 0;
5133	for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5134	int count = columnCounts[iBlock];
5135	columnCounts[iBlock] = kBlock;
5136	kBlock += count;
5137	}
5138	for (iBlock = 0; iBlock < numberElementBlocks; iBlock++) {
5139	int iColumnBlock = blockInfo[iBlock].columnBlock;
5140	whichBlock[columnCounts[iColumnBlock]] = iBlock;
5141	columnCounts[iColumnBlock]++;
5142	}
5143	for (iBlock = numberColumnBlocks - 1; iBlock >= 0; iBlock--)
5144	columnCounts[iBlock + 1] = columnCounts[iBlock];
5145	columnCounts[0] = 0;
5146	for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5147	int count = columnCounts[iBlock + 1] - columnCounts[iBlock];
5148	if (count == 1) {
5149	int kBlock = whichBlock[columnCounts[iBlock]];
5150	int iRowBlock = blockInfo[kBlock].rowBlock;
5151	if (iRowBlock == masterRowBlock) {
5152	if (masterColumnBlock == -1) {
5153	masterColumnBlock = iBlock;
5154	} else {
5155	// Can't decode
5156	masterColumnBlock = -2;
5157	break;
5158	}
5159	}
5160	}
5161	}
5162	if (masterRowBlock < 0 \|\| masterColumnBlock == -2) {
5163	// What now
5164	abort();
5165	}
5166	delete[] rowCounts;
5167	// create all data
5168	const CoinPackedMatrix *top = new const CoinPackedMatrix [numberColumnBlocks];
5169	ClpSimplex *sub = new ClpSimplex[numberBlocks];
5170	ClpSimplex master;
5171	// Set offset
5172	master.setObjectiveOffset(model->objectiveOffset());
5173	bool reducePrint = logLevel() == 7;
5174	if (reducePrint) {
5175	// special
5176	setLogLevel(1);
5177	master.setLogLevel(1);
5178	}
5179	kBlock = 0;
5180	int masterBlock = -1;
5181	for (iBlock = 0; iBlock < numberColumnBlocks; iBlock++) {
5182	top[kBlock] = NULL;
5183	int start = columnCounts[iBlock];
5184	int end = columnCounts[iBlock + 1];
5185	assert(end - start <= 2);
5186	for (int j = start; j < end; j++) {
5187	int jBlock = whichBlock[j];
5188	int iRowBlock = blockInfo[jBlock].rowBlock;
5189	int iColumnBlock = blockInfo[jBlock].columnBlock;
5190	assert(iColumnBlock == iBlock);
5191	if (iColumnBlock != masterColumnBlock && iRowBlock == masterRowBlock) {
5192	// top matrix
5193	top[kBlock] = model->coinBlock(jBlock)->packedMatrix();
5194	} else {
5195	const CoinPackedMatrix *matrix
5196	= model->coinBlock(jBlock)->packedMatrix();
5197	// Get pointers to arrays
5198	const double *rowLower;
5199	const double *rowUpper;
5200	const double *columnLower;
5201	const double *columnUpper;
5202	const double *objective;
5203	model->block(iRowBlock, iColumnBlock, rowLower, rowUpper,
5204	columnLower, columnUpper, objective);
5205	if (iColumnBlock != masterColumnBlock) {
5206	// diagonal block
5207	sub[kBlock].loadProblem(*matrix, columnLower, columnUpper,
5208	objective, rowLower, rowUpper);
5209	if (true) {
5210	double *lower = sub[kBlock].columnLower();
5211	double *upper = sub[kBlock].columnUpper();
5212	int n = sub[kBlock].numberColumns();
5213	for (int i = 0; i < n; i++) {
5214	lower[i] = CoinMax(-1.0e8, lower[i]);
5215	upper[i] = CoinMin(1.0e8, upper[i]);
5216	}
5217	}
5218	if (optimizationDirection_ < 0.0) {
5219	double *obj = sub[kBlock].objective();
5220	int n = sub[kBlock].numberColumns();
5221	for (int i = 0; i < n; i++)
5222	obj[i] = -obj[i];
5223	}
5224	if (this->factorizationFrequency() == 200) {
5225	// User did not touch preset
5226	sub[kBlock].defaultFactorizationFrequency();
5227	} else {
5228	// make sure model has correct value
5229	sub[kBlock].setFactorizationFrequency(this->factorizationFrequency());
5230	}
5231	sub[kBlock].setPerturbation(50);
5232	// Set columnCounts to be diagonal block index for cleanup
5233	columnCounts[kBlock] = jBlock;
5234	} else {
5235	// master
5236	masterBlock = jBlock;
5237	master.loadProblem(*matrix, columnLower, columnUpper,
5238	objective, rowLower, rowUpper);
5239	if (optimizationDirection_ < 0.0) {
5240	double *obj = master.objective();
5241	int n = master.numberColumns();
5242	for (int i = 0; i < n; i++)
5243	obj[i] = -obj[i];
5244	}
5245	}
5246	}
5247	}
5248	if (iBlock != masterColumnBlock)
5249	kBlock++;
5250	}
5251	delete[] whichBlock;
5252	delete[] blockInfo;
5253	// For now master must have been defined (does not have to have columns)
5254	assert(master.numberRows());
5255	assert(masterBlock >= 0);
5256	int numberMasterRows = master.numberRows();
5257	// Overkill in terms of space
5258	int spaceNeeded = CoinMax(numberBlocks * (numberMasterRows + 1),
5259	2 * numberMasterRows);
5260	int *rowAdd = new int[spaceNeeded];
5261	double *elementAdd = new double[spaceNeeded];
5262	spaceNeeded = numberBlocks;
5263	CoinBigIndex *columnAdd = new CoinBigIndex[spaceNeeded + 1];
5264	double *objective = new double[spaceNeeded];
5265	// Add in costed slacks
5266	int firstArtificial = master.numberColumns();
5267	int lastArtificial = firstArtificial;
5268	if (true) {
5269	const double *lower = master.rowLower();
5270	const double *upper = master.rowUpper();
5271	int kCol = 0;
5272	for (int iRow = 0; iRow < numberMasterRows; iRow++) {
5273	if (lower[iRow] > -1.0e10) {
5274	rowAdd[kCol] = iRow;
5275	elementAdd[kCol++] = 1.0;
5276	}
5277	if (upper[iRow] < 1.0e10) {
5278	rowAdd[kCol] = iRow;
5279	elementAdd[kCol++] = -1.0;
5280	}
5281	}
5282	if (kCol > spaceNeeded) {
5283	spaceNeeded = kCol;
5284	delete[] columnAdd;
5285	delete[] objective;
5286	columnAdd = new CoinBigIndex[spaceNeeded + 1];
5287	objective = new double[spaceNeeded];
5288	}
5289	for (int i = 0; i < kCol; i++) {
5290	columnAdd[i] = i;
5291	objective[i] = 1.0e13;
5292	}
5293	columnAdd[kCol] = kCol;
5294	master.addColumns(kCol, NULL, NULL, objective,
5295	columnAdd, rowAdd, elementAdd);
5296	lastArtificial = master.numberColumns();
5297	}
5298	int maxPass = options.independentOption(2);
5299	if (maxPass < 2)
5300	maxPass = 100;
5301	int iPass;
5302	double lastObjective = 1.0e31;
5303	// Create convexity rows for proposals
5304	int numberMasterColumns = master.numberColumns();
5305	master.resize(numberMasterRows + numberBlocks, numberMasterColumns);
5306	if (this->factorizationFrequency() == 200) {
5307	// User did not touch preset
5308	master.defaultFactorizationFrequency();
5309	} else {
5310	// make sure model has correct value
5311	master.setFactorizationFrequency(this->factorizationFrequency());
5312	}
5313	master.setPerturbation(50);
5314	// Arrays to say which block and when created
5315	int maximumColumns = 2 * numberMasterRows + 10 * numberBlocks;
5316	whichBlock = new int[maximumColumns];
5317	int *when = new int[maximumColumns];
5318	int numberColumnsGenerated = numberBlocks;
5319	// fill in rhs and add in artificials
5320	{
5321	double *rowLower = master.rowLower();
5322	double *rowUpper = master.rowUpper();
5323	int iBlock;
5324	columnAdd[0] = 0;
5325	for (iBlock = 0; iBlock < numberBlocks; iBlock++) {
5326	int iRow = iBlock + numberMasterRows;
5327	;
5328	rowLower[iRow] = 1.0;
5329	rowUpper[iRow] = 1.0;
5330	rowAdd[iBlock] = iRow;
5331	elementAdd[iBlock] = 1.0;
5332	objective[iBlock] = 1.0e13;
5333	columnAdd[iBlock + 1] = iBlock + 1;
5334	when[iBlock] = -1;
5335	whichBlock[iBlock] = iBlock;
5336	}
5337	master.addColumns(numberBlocks, NULL, NULL, objective,
5338	columnAdd, rowAdd, elementAdd);
5339	}
5340	char generalPrint[200];
5341	// and resize matrix to double check clp will be happy
5342	//master.matrix()->setDimensions(numberMasterRows+numberBlocks,
5343	// numberMasterColumns+numberBlocks);
5344	sprintf(generalPrint, "Time to decompose %.2f seconds", CoinCpuTime() - time1);
5345	handler_->message(CLP_GENERAL, messages_)
5346	<< generalPrint
5347	<< CoinMessageEol;
5348	#ifdef ABC_INHERIT
5349	//AbcSimplex abcMaster;
5350	//if (!this->abcState())
5351	//setAbcState(1);
5352	int numberCpu = CoinMin((this->abcState() & 15), 4);
5353	CoinPthreadStuff threadInfo(numberCpu, clp_parallelManager);
5354	master.setAbcState(this->abcState());
5355	//AbcSimplex * tempMaster=master.dealWithAbc(2,10,true);
5356	//abcMaster=*tempMaster;
5357	//delete tempMaster;
5358	//abcMaster.startThreads(numberCpu);
5359	//#define master abcMaster
5360	#endif
5361	for (iPass = 0; iPass < maxPass; iPass++) {
5362	sprintf(generalPrint, "Start of pass %d", iPass);
5363	handler_->message(CLP_GENERAL, messages_)
5364	<< generalPrint
5365	<< CoinMessageEol;
5366	// Solve master - may be infeasible
5367	//master.scaling(0);
5368	if (0) {
5369	master.writeMps("yy.mps");
5370	}
5371	// Correct artificials
5372	double sumArtificials = 0.0;
5373	if (iPass) {
5374	double *upper = master.columnUpper();
5375	double *solution = master.primalColumnSolution();
5376	double *obj = master.objective();
5377	sumArtificials = 0.0;
5378	for (int i = firstArtificial; i < lastArtificial; i++) {
5379	sumArtificials += solution[i];
5380	//assert (solution[i]>-1.0e-2);
5381	if (solution[i] < 1.0e-6) {
5382	#if 0
5383	// Could take out
5384	obj[i] = 0.0;
5385	upper[i] = 0.0;
5386	#else
5387	obj[i] = 1.0e7;
5388	upper[i] = 1.0e-1;
5389	#endif
5390	solution[i] = 0.0;
5391	master.setColumnStatus(i, isFixed);
5392	} else {
5393	upper[i] = solution[i] + 1.0e-5 * (1.0 + solution[i]);
5394	}
5395	}
5396	sprintf(generalPrint, "Sum of artificials before solve is %g", sumArtificials);
5397	handler_->message(CLP_GENERAL, messages_)
5398	<< generalPrint
5399	<< CoinMessageEol;
5400	}
5401	// scale objective to be reasonable
5402	double scaleFactor = master.scaleObjective(-1.0e9);
5403	{
5404	double *dual = master.dualRowSolution();
5405	int n = master.numberRows();
5406	memset(dual, 0, n * sizeof(double));
5407	double *solution = master.primalColumnSolution();
5408	master.clpMatrix()->times(1.0, solution, dual);
5409	double sum = 0.0;
5410	double *lower = master.rowLower();
5411	double *upper = master.rowUpper();
5412	for (int iRow = 0; iRow < n; iRow++) {
5413	double value = dual[iRow];
5414	if (value > upper[iRow])
5415	sum += value - upper[iRow];
5416	else if (value < lower[iRow])
5417	sum -= value - lower[iRow];
5418	}
5419	printf("** suminf %g\n", sum);
5420	lower = master.columnLower();
5421	upper = master.columnUpper();
5422	n = master.numberColumns();
5423	for (int iColumn = 0; iColumn < n; iColumn++) {
5424	double value = solution[iColumn];
5425	if (value > upper[iColumn] + 1.0e-5)
5426	sum += value - upper[iColumn];
5427	else if (value < lower[iColumn] - 1.0e-5)
5428	sum -= value - lower[iColumn];
5429	}
5430	printf("** suminf %g\n", sum);
5431	}
5432	#ifdef ABC_INHERIT
5433	master.dealWithAbc(1, 0, true);
5434	#else
5435	master.primal();
5436	#endif
5437	//master.primal(1);
5438	// Correct artificials
5439	sumArtificials = 0.0;
5440	{
5441	double *solution = master.primalColumnSolution();
5442	for (int i = firstArtificial; i < lastArtificial; i++) {
5443	sumArtificials += solution[i];
5444	}
5445	printf("** Sum of artificials after solve is %g\n", sumArtificials);
5446	}
5447	master.scaleObjective(scaleFactor);
5448	int problemStatus = master.status(); // do here as can change (delcols)
5449	if (problemStatus == 2 && master.numberColumns()) {
5450	#ifdef ABC_INHERIT
5451	master.dealWithAbc(1, 1, true);
5452	#else
5453	master.primal(1);
5454	#endif
5455	//master.primal(1);
5456	if (problemStatus == 2)