Context Navigation

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

Last change on this file since 1610 was 1610, checked in by forrest, 11 years ago
changes for more flexibility in ClpSolve?
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 201.6 KB

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

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

Download in other formats: