Context Navigation

source: trunk/Clp/src/ClpSolve.cpp @ 1525

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago
Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 200.7 KB

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

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

Download in other formats: