Context Navigation

ClpSimplexDual.cpp @ 1780

Last change on this file since 1780 was 1780, checked in by forrest, 9 years ago
stuff for parametrics
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 325.6 KB

Line
1	/* $Id: ClpSimplexDual.cpp 1780 2011-08-11 07:15:40Z forrest $ */
2	// Copyright (C) 2002, International Business Machines
3	// Corporation and others. All Rights Reserved.
4	// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7	/* Notes on implementation of dual simplex algorithm.
8
9	When dual feasible:
10
11	If primal feasible, we are optimal. Otherwise choose an infeasible
12	basic variable to leave basis (normally going to nearest bound) (B). We
13	now need to find an incoming variable which will leave problem
14	dual feasible so we get the row of the tableau corresponding to
15	the basic variable (with the correct sign depending if basic variable
16	above or below feasibility region - as that affects whether reduced
17	cost on outgoing variable has to be positive or negative).
18
19	We now perform a ratio test to determine which incoming variable will
20	preserve dual feasibility (C). If no variable found then problem
21	is infeasible (in primal sense). If there is a variable, we then
22	perform pivot and repeat. Trivial?
23
24	-------------------------------------------
25
26	A) How do we get dual feasible? If all variables have bounds then
27	it is trivial to get feasible by putting non-basic variables to
28	correct bounds. OSL did not have a phase 1/phase 2 approach but
29	instead effectively put fake bounds on variables and this is the
30	approach here, although I had hoped to make it cleaner.
31
32	If there is a weight of X on getting dual feasible:
33	Non-basic variables with negative reduced costs are put to
34	lesser of their upper bound and their lower bound + X.
35	Similarly, mutatis mutandis, for positive reduced costs.
36
37	Free variables should normally be in basis, otherwise I have
38	coding which may be able to come out (and may not be correct).
39
40	In OSL, this weight was changed heuristically, here at present
41	it is only increased if problem looks finished. If problem is
42	feasible I check for unboundedness. If not unbounded we
43	could play with going into primal. As long as weights increase
44	any algorithm would be finite.
45
46	B) Which outgoing variable to choose is a virtual base class.
47	For difficult problems steepest edge is preferred while for
48	very easy (large) problems we will need partial scan.
49
50	C) Sounds easy, but this is hardest part of algorithm.
51	1) Instead of stopping at first choice, we may be able
52	to flip that variable to other bound and if objective
53	still improving choose again. These mini iterations can
54	increase speed by orders of magnitude but we may need to
55	go to more of a bucket choice of variable rather than looking
56	at them one by one (for speed).
57	2) Accuracy. Reduced costs may be of wrong sign but less than
58	tolerance. Pivoting on these makes objective go backwards.
59	OSL modified cost so a zero move was made, Gill et al
60	(in primal analogue) modified so a strictly positive move was
61	made. It is not quite as neat in dual but that is what we
62	try and do. The two problems are that re-factorizations can
63	change reduced costs above and below tolerances and that when
64	finished we need to reset costs and try again.
65	3) Degeneracy. Gill et al helps but may not be enough. We
66	may need more. Also it can improve speed a lot if we perturb
67	the costs significantly.
68
69	References:
70	Forrest and Goldfarb, Steepest-edge simplex algorithms for
71	linear programming - Mathematical Programming 1992
72	Forrest and Tomlin, Implementing the simplex method for
73	the Optimization Subroutine Library - IBM Systems Journal 1992
74	Gill, Murray, Saunders, Wright A Practical Anti-Cycling
75	Procedure for Linear and Nonlinear Programming SOL report 1988
76
77
78	TODO:
79
80	a) Better recovery procedures. At present I never check on forward
81	progress. There is checkpoint/restart with reducing
82	re-factorization frequency, but this is only on singular
83	factorizations.
84	b) Fast methods for large easy problems (and also the option for
85	the code to automatically choose which method).
86	c) We need to be able to stop in various ways for OSI - this
87	is fairly easy.
88
89	*/
90	#ifdef COIN_DEVELOP
91	#undef COIN_DEVELOP
92	#define COIN_DEVELOP 2
93	#endif
94
95	#include "CoinPragma.hpp"
96
97	#include <math.h>
98
99	#include "CoinHelperFunctions.hpp"
100	#include "ClpHelperFunctions.hpp"
101	#include "ClpSimplexDual.hpp"
102	#include "ClpEventHandler.hpp"
103	#include "ClpFactorization.hpp"
104	#include "CoinPackedMatrix.hpp"
105	#include "CoinIndexedVector.hpp"
106	#include "CoinFloatEqual.hpp"
107	#include "ClpDualRowDantzig.hpp"
108	#include "ClpMessage.hpp"
109	#include "ClpLinearObjective.hpp"
110	#include <cfloat>
111	#include <cassert>
112	#include <string>
113	#include <stdio.h>
114	#include <iostream>
115	//#define CLP_DEBUG 1
116	// To force to follow another run put logfile name here and define
117	//#define FORCE_FOLLOW
118	#ifdef FORCE_FOLLOW
119	static FILE * fpFollow = NULL;
120	static char * forceFile = "old.log";
121	static int force_in = -1;
122	static int force_out = -1;
123	static int force_iteration = 0;
124	#endif
125	//#define VUB
126	#ifdef VUB
127	extern int * vub;
128	extern int * toVub;
129	extern int * nextDescendent;
130	#endif
131	#ifdef NDEBUG
132	#define NDEBUG_CLP
133	#endif
134	#ifndef CLP_INVESTIGATE
135	#define NDEBUG_CLP
136	#endif
137	// dual
138
139	/* *** Method
140	This is a vanilla version of dual simplex.
141
142	It tries to be a single phase approach with a weight of 1.0 being
143	given to getting optimal and a weight of dualBound_ being
144	given to getting dual feasible. In this version I have used the
145	idea that this weight can be thought of as a fake bound. If the
146	distance between the lower and upper bounds on a variable is less
147	than the feasibility weight then we are always better off flipping
148	to other bound to make dual feasible. If the distance is greater
149	then we make up a fake bound dualBound_ away from one bound.
150	If we end up optimal or primal infeasible, we check to see if
151	bounds okay. If so we have finished, if not we increase dualBound_
152	and continue (after checking if unbounded). I am undecided about
153	free variables - there is coding but I am not sure about it. At
154	present I put them in basis anyway.
155
156	The code is designed to take advantage of sparsity so arrays are
157	seldom zeroed out from scratch or gone over in their entirety.
158	The only exception is a full scan to find outgoing variable. This
159	will be changed to keep an updated list of infeasibilities (or squares
160	if steepest edge). Also on easy problems we don't need full scan - just
161	pick first reasonable.
162
163	One problem is how to tackle degeneracy and accuracy. At present
164	I am using the modification of costs which I put in OSL and which was
165	extended by Gill et al. I am still not sure of the exact details.
166
167	The flow of dual is three while loops as follows:
168
169	while (not finished) {
170
171	while (not clean solution) {
172
173	Factorize and/or clean up solution by flipping variables so
174	dual feasible. If looks finished check fake dual bounds.
175	Repeat until status is iterating (-1) or finished (0,1,2)
176
177	}
178
179	while (status==-1) {
180
181	Iterate until no pivot in or out or time to re-factorize.
182
183	Flow is:
184
185	choose pivot row (outgoing variable). if none then
186	we are primal feasible so looks as if done but we need to
187	break and check bounds etc.
188
189	Get pivot row in tableau
190
191	Choose incoming column. If we don't find one then we look
192	primal infeasible so break and check bounds etc. (Also the
193	pivot tolerance is larger after any iterations so that may be
194	reason)
195
196	If we do find incoming column, we may have to adjust costs to
197	keep going forwards (anti-degeneracy). Check pivot will be stable
198	and if unstable throw away iteration (we will need to implement
199	flagging of basic variables sometime) and break to re-factorize.
200	If minor error re-factorize after iteration.
201
202	Update everything (this may involve flipping variables to stay
203	dual feasible.
204
205	}
206
207	}
208
209	At present we never check we are going forwards. I overdid that in
210	OSL so will try and make a last resort.
211
212	Needs partial scan pivot out option.
213	Needs dantzig, uninitialized and full steepest edge options (can still
214	use partial scan)
215
216	May need other anti-degeneracy measures, especially if we try and use
217	loose tolerances as a way to solve in fewer iterations.
218
219	I like idea of dynamic scaling. This gives opportunity to decouple
220	different implications of scaling for accuracy, iteration count and
221	feasibility tolerance.
222
223	*/
224	#define CLEAN_FIXED 0
225	// Startup part of dual (may be extended to other algorithms)
226	int
227	ClpSimplexDual::startupSolve(int ifValuesPass, double * saveDuals, int startFinishOptions)
228	{
229	// If values pass then save given duals round check solution
230	// sanity check
231	// initialize - no values pass and algorithm_ is -1
232	// put in standard form (and make row copy)
233	// create modifiable copies of model rim and do optional scaling
234	// If problem looks okay
235	// Do initial factorization
236	// If user asked for perturbation - do it
237	numberFake_ = 0; // Number of variables at fake bounds
238	numberChanged_ = 0; // Number of variables with changed costs
239	if (!startup(0, startFinishOptions)) {
240	int usePrimal = 0;
241	// looks okay
242	// Superbasic variables not allowed
243	// If values pass then scale pi
244	if (ifValuesPass) {
245	if (problemStatus_ && perturbation_ < 100)
246	usePrimal = perturb();
247	int i;
248	if (scalingFlag_ > 0) {
249	for (i = 0; i < numberRows_; i++) {
250	dual_[i] = saveDuals[i] * inverseRowScale_[i];
251	}
252	} else {
253	CoinMemcpyN(saveDuals, numberRows_, dual_);
254	}
255	// now create my duals
256	for (i = 0; i < numberRows_; i++) {
257	// slack
258	double value = dual_[i];
259	value += rowObjectiveWork_[i];
260	saveDuals[i+numberColumns_] = value;
261	}
262	CoinMemcpyN(objectiveWork_, numberColumns_, saveDuals);
263	transposeTimes(-1.0, dual_, saveDuals);
264	// make reduced costs okay
265	for (i = 0; i < numberColumns_; i++) {
266	if (getStatus(i) == atLowerBound) {
267	if (saveDuals[i] < 0.0) {
268	//if (saveDuals[i]<-1.0e-3)
269	//printf("bad dj at lb %d %g\n",i,saveDuals[i]);
270	saveDuals[i] = 0.0;
271	}
272	} else if (getStatus(i) == atUpperBound) {
273	if (saveDuals[i] > 0.0) {
274	//if (saveDuals[i]>1.0e-3)
275	//printf("bad dj at ub %d %g\n",i,saveDuals[i]);
276	saveDuals[i] = 0.0;
277	}
278	}
279	}
280	CoinMemcpyN(saveDuals, (numberColumns_ + numberRows_), dj_);
281	// set up possible ones
282	for (i = 0; i < numberRows_ + numberColumns_; i++)
283	clearPivoted(i);
284	int iRow;
285	for (iRow = 0; iRow < numberRows_; iRow++) {
286	int iPivot = pivotVariable_[iRow];
287	if (fabs(saveDuals[iPivot]) > dualTolerance_) {
288	if (getStatus(iPivot) != isFree)
289	setPivoted(iPivot);
290	}
291	}
292	} else if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
293	// set up possible ones
294	for (int i = 0; i < numberRows_ + numberColumns_; i++)
295	clearPivoted(i);
296	int iRow;
297	for (iRow = 0; iRow < numberRows_; iRow++) {
298	int iPivot = pivotVariable_[iRow];
299	if (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]) {
300	setPivoted(iPivot);
301	}
302	}
303	}
304
305	double objectiveChange;
306	assert (!numberFake_);
307	assert (numberChanged_ == 0);
308	if (!numberFake_) // if nonzero then adjust
309	changeBounds(1, NULL, objectiveChange);
310
311	if (!ifValuesPass) {
312	// Check optimal
313	if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_)
314	problemStatus_ = 0;
315	}
316	if (problemStatus_ < 0 && perturbation_ < 100) {
317	bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
318	if (!inCbcOrOther)
319	usePrimal = perturb();
320	// Can't get here if values pass
321	gutsOfSolution(NULL, NULL);
322	#ifdef CLP_INVESTIGATE
323	if (numberDualInfeasibilities_)
324	printf("ZZZ %d primal %d dual - sumdinf %g\n",
325	numberPrimalInfeasibilities_,
326	numberDualInfeasibilities_, sumDualInfeasibilities_);
327	#endif
328	if (handler_->logLevel() > 2) {
329	handler_->message(CLP_SIMPLEX_STATUS, messages_)
330	<< numberIterations_ << objectiveValue();
331	handler_->printing(sumPrimalInfeasibilities_ > 0.0)
332	<< sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
333	handler_->printing(sumDualInfeasibilities_ > 0.0)
334	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
335	handler_->printing(numberDualInfeasibilitiesWithoutFree_
336	< numberDualInfeasibilities_)
337	<< numberDualInfeasibilitiesWithoutFree_;
338	handler_->message() << CoinMessageEol;
339	}
340	if (inCbcOrOther) {
341	if (numberPrimalInfeasibilities_) {
342	usePrimal = perturb();
343	if (perturbation_ >= 101) {
344	computeDuals(NULL);
345	//gutsOfSolution(NULL,NULL);
346	checkDualSolution(); // recompute objective
347	}
348	} else if (numberDualInfeasibilities_) {
349	problemStatus_ = 10;
350	if ((moreSpecialOptions_ & 32) != 0 && false)
351	problemStatus_ = 0; // say optimal!!
352	#if COIN_DEVELOP>2
353
354	printf("returning at %d\n", __LINE__);
355	#endif
356	return 1; // to primal
357	}
358	}
359	} else if (!ifValuesPass) {
360	gutsOfSolution(NULL, NULL);
361	// double check
362	if (numberDualInfeasibilities_ \|\| numberPrimalInfeasibilities_)
363	problemStatus_ = -1;
364	}
365	if (usePrimal) {
366	problemStatus_ = 10;
367	#if COIN_DEVELOP>2
368	printf("returning to use primal (no obj) at %d\n", __LINE__);
369	#endif
370	}
371	return usePrimal;
372	} else {
373	return 1;
374	}
375	}
376	void
377	ClpSimplexDual::finishSolve(int startFinishOptions)
378	{
379	assert (problemStatus_ \|\| !sumPrimalInfeasibilities_);
380
381	// clean up
382	finish(startFinishOptions);
383	}
384	//#define CLP_REPORT_PROGRESS
385	#ifdef CLP_REPORT_PROGRESS
386	static int ixxxxxx = 0;
387	static int ixxyyyy = 90;
388	#endif
389	#ifdef CLP_INVESTIGATE_SERIAL
390	static int z_reason[7] = {0, 0, 0, 0, 0, 0, 0};
391	static int z_thinks = -1;
392	#endif
393	void
394	ClpSimplexDual::gutsOfDual(int ifValuesPass, double * & saveDuals, int initialStatus,
395	ClpDataSave & data)
396	{
397	#ifdef CLP_INVESTIGATE_SERIAL
398	z_reason[0]++;
399	z_thinks = -1;
400	int nPivots = 9999;
401	#endif
402	double largestPrimalError = 0.0;
403	double largestDualError = 0.0;
404	// Start can skip some things in transposeTimes
405	specialOptions_ \|= 131072;
406	int lastCleaned = 0; // last time objective or bounds cleaned up
407
408	// This says whether to restore things etc
409	// startup will have factorized so can skip
410	int factorType = 0;
411	// Start check for cycles
412	progress_.startCheck();
413	// Say change made on first iteration
414	changeMade_ = 1;
415	// Say last objective infinite
416	//lastObjectiveValue_=-COIN_DBL_MAX;
417	progressFlag_ = 0;
418	/*
419	Status of problem:
420	0 - optimal
421	1 - infeasible
422	2 - unbounded
423	-1 - iterating
424	-2 - factorization wanted
425	-3 - redo checking without factorization
426	-4 - looks infeasible
427	*/
428	while (problemStatus_ < 0) {
429	int iRow, iColumn;
430	// clear
431	for (iRow = 0; iRow < 4; iRow++) {
432	rowArray_[iRow]->clear();
433	}
434
435	for (iColumn = 0; iColumn < 2; iColumn++) {
436	columnArray_[iColumn]->clear();
437	}
438
439	// give matrix (and model costs and bounds a chance to be
440	// refreshed (normally null)
441	matrix_->refresh(this);
442	// If getting nowhere - why not give it a kick
443	// does not seem to work too well - do some more work
444	if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_)
445	&& initialStatus != 10) {
446	perturb();
447	// Can't get here if values pass
448	gutsOfSolution(NULL, NULL);
449	if (handler_->logLevel() > 2) {
450	handler_->message(CLP_SIMPLEX_STATUS, messages_)
451	<< numberIterations_ << objectiveValue();
452	handler_->printing(sumPrimalInfeasibilities_ > 0.0)
453	<< sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
454	handler_->printing(sumDualInfeasibilities_ > 0.0)
455	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
456	handler_->printing(numberDualInfeasibilitiesWithoutFree_
457	< numberDualInfeasibilities_)
458	<< numberDualInfeasibilitiesWithoutFree_;
459	handler_->message() << CoinMessageEol;
460	}
461	}
462	// see if in Cbc etc
463	bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
464	#if 0
465	bool gotoPrimal = false;
466	if (inCbcOrOther && numberIterations_ > disasterArea_ + numberRows_ &&
467	numberDualInfeasibilitiesWithoutFree_ && largestDualError_ > 1.0e-1) {
468	if (!disasterArea_) {
469	printf("trying all slack\n");
470	// try all slack basis
471	allSlackBasis(true);
472	disasterArea_ = 2 * numberRows_;
473	} else {
474	printf("going to primal\n");
475	// go to primal
476	gotoPrimal = true;
477	allSlackBasis(true);
478	}
479	}
480	#endif
481	bool disaster = false;
482	if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
483	disasterArea_->saveInfo();
484	disaster = true;
485	}
486	// may factorize, checks if problem finished
487	statusOfProblemInDual(lastCleaned, factorType, saveDuals, data,
488	ifValuesPass);
489	largestPrimalError = CoinMax(largestPrimalError, largestPrimalError_);
490	largestDualError = CoinMax(largestDualError, largestDualError_);
491	if (disaster)
492	problemStatus_ = 3;
493	// If values pass then do easy ones on first time
494	if (ifValuesPass &&
495	progress_.lastIterationNumber(0) < 0 && saveDuals) {
496	doEasyOnesInValuesPass(saveDuals);
497	}
498
499	// Say good factorization
500	factorType = 1;
501	if (data.sparseThreshold_) {
502	// use default at present
503	factorization_->sparseThreshold(0);
504	factorization_->goSparse();
505	}
506
507	// exit if victory declared
508	if (problemStatus_ >= 0)
509	break;
510
511	// test for maximum iterations
512	if (hitMaximumIterations() \|\| (ifValuesPass == 2 && !saveDuals)) {
513	problemStatus_ = 3;
514	break;
515	}
516	if (ifValuesPass && !saveDuals) {
517	// end of values pass
518	ifValuesPass = 0;
519	int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
520	if (status >= 0) {
521	problemStatus_ = 5;
522	secondaryStatus_ = ClpEventHandler::endOfValuesPass;
523	break;
524	}
525	}
526	// Check event
527	{
528	int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
529	if (status >= 0) {
530	problemStatus_ = 5;
531	secondaryStatus_ = ClpEventHandler::endOfFactorization;
532	break;
533	}
534	}
535	// Do iterations
536	int returnCode = whileIterating(saveDuals, ifValuesPass);
537	if (problemStatus_ == 1 && (progressFlag_&8) != 0 &&
538	fabs(objectiveValue_) > 1.0e10 )
539	problemStatus_ = 10; // infeasible - but has looked feasible
540	#ifdef CLP_INVESTIGATE_SERIAL
541	nPivots = factorization_->pivots();
542	#endif
543	if (!problemStatus_ && factorization_->pivots())
544	computeDuals(NULL); // need to compute duals
545	if (returnCode == -2)
546	factorType = 3;
547	}
548	#ifdef CLP_INVESTIGATE_SERIAL
549	// NOTE - can fail if parallel
550	if (z_thinks != -1) {
551	assert (z_thinks < 4);
552	if ((!factorization_->pivots() && nPivots < 20) && z_thinks >= 0 && z_thinks < 2)
553	z_thinks += 4;
554	z_reason[1+z_thinks]++;
555	}
556	if ((z_reason[0] % 1000) == 0) {
557	printf("Reason");
558	for (int i = 0; i < 7; i++)
559	printf(" %d", z_reason[i]);
560	printf("\n");
561	}
562	#endif
563	// Stop can skip some things in transposeTimes
564	specialOptions_ &= ~131072;
565	largestPrimalError_ = largestPrimalError;
566	largestDualError_ = largestDualError;
567	}
568	int
569	ClpSimplexDual::dual(int ifValuesPass, int startFinishOptions)
570	{
571	bestObjectiveValue_ = -COIN_DBL_MAX;
572	algorithm_ = -1;
573	moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
574	// save data
575	ClpDataSave data = saveData();
576	double * saveDuals = NULL;
577	int saveDont = dontFactorizePivots_;
578	if ((specialOptions_ & 2048) == 0)
579	dontFactorizePivots_ = 0;
580	else if(!dontFactorizePivots_)
581	dontFactorizePivots_ = 20;
582	if (ifValuesPass) {
583	saveDuals = new double [numberRows_+numberColumns_];
584	CoinMemcpyN(dual_, numberRows_, saveDuals);
585	}
586	if (alphaAccuracy_ != -1.0)
587	alphaAccuracy_ = 1.0;
588	int returnCode = startupSolve(ifValuesPass, saveDuals, startFinishOptions);
589	// Save so can see if doing after primal
590	int initialStatus = problemStatus_;
591	if (!returnCode && !numberDualInfeasibilities_ &&
592	!numberPrimalInfeasibilities_ && perturbation_ < 101) {
593	returnCode = 1; // to skip gutsOfDual
594	problemStatus_ = 0;
595	}
596
597	if (!returnCode)
598	gutsOfDual(ifValuesPass, saveDuals, initialStatus, data);
599	if (!problemStatus_) {
600	// see if cutoff reached
601	double limit = 0.0;
602	getDblParam(ClpDualObjectiveLimit, limit);
603	if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
604	limit + 1.0e-7 + 1.0e-8 * fabs(limit)) {
605	// actually infeasible on objective
606	problemStatus_ = 1;
607	secondaryStatus_ = 1;
608	}
609	}
610	// If infeasible but primal errors - try dual
611	if (problemStatus_==1 && numberPrimalInfeasibilities_) {
612	bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
613	double factor = (!inCbcOrOther) ? 1.0 : 0.3;
614	double averageInfeasibility = sumPrimalInfeasibilities_/
615	static_cast<double>(numberPrimalInfeasibilities_);
616	if (averageInfeasibility<factor*largestPrimalError_)
617	problemStatus_= 10;
618	}
619
620	if (problemStatus_ == 10)
621	startFinishOptions \|= 1;
622	finishSolve(startFinishOptions);
623	delete [] saveDuals;
624
625	// Restore any saved stuff
626	restoreData(data);
627	dontFactorizePivots_ = saveDont;
628	if (problemStatus_ == 3)
629	objectiveValue_ = CoinMax(bestObjectiveValue_, objectiveValue_ - bestPossibleImprovement_);
630	return problemStatus_;
631	}
632	// old way
633	#if 0
634	int ClpSimplexDual::dual (int ifValuesPass , int startFinishOptions)
635	{
636	algorithm_ = -1;
637
638	// save data
639	ClpDataSave data = saveData();
640	// Save so can see if doing after primal
641	int initialStatus = problemStatus_;
642
643	// If values pass then save given duals round check solution
644	double * saveDuals = NULL;
645	if (ifValuesPass) {
646	saveDuals = new double [numberRows_+numberColumns_];
647	CoinMemcpyN(dual_, numberRows_, saveDuals);
648	}
649	// sanity check
650	// initialize - no values pass and algorithm_ is -1
651	// put in standard form (and make row copy)
652	// create modifiable copies of model rim and do optional scaling
653	// If problem looks okay
654	// Do initial factorization
655	// If user asked for perturbation - do it
656	if (!startup(0, startFinishOptions)) {
657	// looks okay
658	// Superbasic variables not allowed
659	// If values pass then scale pi
660	if (ifValuesPass) {
661	if (problemStatus_ && perturbation_ < 100)
662	perturb();
663	int i;
664	if (scalingFlag_ > 0) {
665	for (i = 0; i < numberRows_; i++) {
666	dual_[i] = saveDuals[i] * inverseRowScale_[i];
667	}
668	} else {
669	CoinMemcpyN(saveDuals, numberRows_, dual_);
670	}
671	// now create my duals
672	for (i = 0; i < numberRows_; i++) {
673	// slack
674	double value = dual_[i];
675	value += rowObjectiveWork_[i];
676	saveDuals[i+numberColumns_] = value;
677	}
678	CoinMemcpyN(objectiveWork_, numberColumns_, saveDuals);
679	transposeTimes(-1.0, dual_, saveDuals);
680	// make reduced costs okay
681	for (i = 0; i < numberColumns_; i++) {
682	if (getStatus(i) == atLowerBound) {
683	if (saveDuals[i] < 0.0) {
684	//if (saveDuals[i]<-1.0e-3)
685	//printf("bad dj at lb %d %g\n",i,saveDuals[i]);
686	saveDuals[i] = 0.0;
687	}
688	} else if (getStatus(i) == atUpperBound) {
689	if (saveDuals[i] > 0.0) {
690	//if (saveDuals[i]>1.0e-3)
691	//printf("bad dj at ub %d %g\n",i,saveDuals[i]);
692	saveDuals[i] = 0.0;
693	}
694	}
695	}
696	CoinMemcpyN(saveDuals, numberColumns_ + numberRows_, dj_);
697	// set up possible ones
698	for (i = 0; i < numberRows_ + numberColumns_; i++)
699	clearPivoted(i);
700	int iRow;
701	for (iRow = 0; iRow < numberRows_; iRow++) {
702	int iPivot = pivotVariable_[iRow];
703	if (fabs(saveDuals[iPivot]) > dualTolerance_) {
704	if (getStatus(iPivot) != isFree)
705	setPivoted(iPivot);
706	}
707	}
708	} else if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
709	// set up possible ones
710	for (int i = 0; i < numberRows_ + numberColumns_; i++)
711	clearPivoted(i);
712	int iRow;
713	for (iRow = 0; iRow < numberRows_; iRow++) {
714	int iPivot = pivotVariable_[iRow];
715	if (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]) {
716	setPivoted(iPivot);
717	}
718	}
719	}
720
721	double objectiveChange;
722	numberFake_ = 0; // Number of variables at fake bounds
723	numberChanged_ = 0; // Number of variables with changed costs
724	changeBounds(1, NULL, objectiveChange);
725
726	int lastCleaned = 0; // last time objective or bounds cleaned up
727
728	if (!ifValuesPass) {
729	// Check optimal
730	if (!numberDualInfeasibilities_ && !numberPrimalInfeasibilities_)
731	problemStatus_ = 0;
732	}
733	if (problemStatus_ < 0 && perturbation_ < 100) {
734	perturb();
735	// Can't get here if values pass
736	gutsOfSolution(NULL, NULL);
737	if (handler_->logLevel() > 2) {
738	handler_->message(CLP_SIMPLEX_STATUS, messages_)
739	<< numberIterations_ << objectiveValue();
740	handler_->printing(sumPrimalInfeasibilities_ > 0.0)
741	<< sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
742	handler_->printing(sumDualInfeasibilities_ > 0.0)
743	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
744	handler_->printing(numberDualInfeasibilitiesWithoutFree_
745	< numberDualInfeasibilities_)
746	<< numberDualInfeasibilitiesWithoutFree_;
747	handler_->message() << CoinMessageEol;
748	}
749	}
750
751	// This says whether to restore things etc
752	// startup will have factorized so can skip
753	int factorType = 0;
754	// Start check for cycles
755	progress_.startCheck();
756	// Say change made on first iteration
757	changeMade_ = 1;
758	/*
759	Status of problem:
760	0 - optimal
761	1 - infeasible
762	2 - unbounded
763	-1 - iterating
764	-2 - factorization wanted
765	-3 - redo checking without factorization
766	-4 - looks infeasible
767	*/
768	while (problemStatus_ < 0) {
769	int iRow, iColumn;
770	// clear
771	for (iRow = 0; iRow < 4; iRow++) {
772	rowArray_[iRow]->clear();
773	}
774
775	for (iColumn = 0; iColumn < 2; iColumn++) {
776	columnArray_[iColumn]->clear();
777	}
778
779	// give matrix (and model costs and bounds a chance to be
780	// refreshed (normally null)
781	matrix_->refresh(this);
782	// If getting nowhere - why not give it a kick
783	// does not seem to work too well - do some more work
784	if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_)
785	&& initialStatus != 10) {
786	perturb();
787	// Can't get here if values pass
788	gutsOfSolution(NULL, NULL);
789	if (handler_->logLevel() > 2) {
790	handler_->message(CLP_SIMPLEX_STATUS, messages_)
791	<< numberIterations_ << objectiveValue();
792	handler_->printing(sumPrimalInfeasibilities_ > 0.0)
793	<< sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
794	handler_->printing(sumDualInfeasibilities_ > 0.0)
795	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
796	handler_->printing(numberDualInfeasibilitiesWithoutFree_
797	< numberDualInfeasibilities_)
798	<< numberDualInfeasibilitiesWithoutFree_;
799	handler_->message() << CoinMessageEol;
800	}
801	}
802	// may factorize, checks if problem finished
803	statusOfProblemInDual(lastCleaned, factorType, saveDuals, data,
804	ifValuesPass);
805	// If values pass then do easy ones on first time
806	if (ifValuesPass &&
807	progress_.lastIterationNumber(0) < 0 && saveDuals) {
808	doEasyOnesInValuesPass(saveDuals);
809	}
810
811	// Say good factorization
812	factorType = 1;
813	if (data.sparseThreshold_) {
814	// use default at present
815	factorization_->sparseThreshold(0);
816	factorization_->goSparse();
817	}
818
819	// exit if victory declared
820	if (problemStatus_ >= 0)
821	break;
822
823	// test for maximum iterations
824	if (hitMaximumIterations() \|\| (ifValuesPass == 2 && !saveDuals)) {
825	problemStatus_ = 3;
826	break;
827	}
828	if (ifValuesPass && !saveDuals) {
829	// end of values pass
830	ifValuesPass = 0;
831	int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
832	if (status >= 0) {
833	problemStatus_ = 5;
834	secondaryStatus_ = ClpEventHandler::endOfValuesPass;
835	break;
836	}
837	}
838	// Check event
839	{
840	int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
841	if (status >= 0) {
842	problemStatus_ = 5;
843	secondaryStatus_ = ClpEventHandler::endOfFactorization;
844	break;
845	}
846	}
847	// Do iterations
848	whileIterating(saveDuals, ifValuesPass);
849	}
850	}
851
852	assert (problemStatus_ \|\| !sumPrimalInfeasibilities_);
853
854	// clean up
855	finish(startFinishOptions);
856	delete [] saveDuals;
857
858	// Restore any saved stuff
859	restoreData(data);
860	return problemStatus_;
861	}
862	#endif
863	//#define CHECK_ACCURACY
864	#ifdef CHECK_ACCURACY
865	static double zzzzzz[100000];
866	#endif
867	/* Reasons to come out:
868	-1 iterations etc
869	-2 inaccuracy
870	-3 slight inaccuracy (and done iterations)
871	+0 looks optimal (might be unbounded - but we will investigate)
872	+1 looks infeasible
873	+3 max iterations
874	*/
875	int
876	ClpSimplexDual::whileIterating(double * & givenDuals, int ifValuesPass)
877	{
878	#ifdef CLP_INVESTIGATE_SERIAL
879	z_thinks = -1;
880	#endif
881	#ifdef CLP_DEBUG
882	int debugIteration = -1;
883	#endif
884	{
885	int i;
886	for (i = 0; i < 4; i++) {
887	rowArray_[i]->clear();
888	}
889	for (i = 0; i < 2; i++) {
890	columnArray_[i]->clear();
891	}
892	}
893	#ifdef CLP_REPORT_PROGRESS
894	double * savePSol = new double [numberRows_+numberColumns_];
895	double * saveDj = new double [numberRows_+numberColumns_];
896	double * saveCost = new double [numberRows_+numberColumns_];
897	unsigned char * saveStat = new unsigned char [numberRows_+numberColumns_];
898	#endif
899	// if can't trust much and long way from optimal then relax
900	if (largestPrimalError_ > 10.0)
901	factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
902	else
903	factorization_->relaxAccuracyCheck(1.0);
904	// status stays at -1 while iterating, >=0 finished, -2 to invert
905	// status -3 to go to top without an invert
906	int returnCode = -1;
907	double saveSumDual = sumDualInfeasibilities_; // so we know to be careful
908
909	#if 0
910	// compute average infeasibility for backward test
911	double averagePrimalInfeasibility = sumPrimalInfeasibilities_ /
912	((double ) (numberPrimalInfeasibilities_ + 1));
913	#endif
914
915	// Get dubious weights
916	CoinBigIndex * dubiousWeights = NULL;
917	#ifdef DUBIOUS_WEIGHTS
918	factorization_->getWeights(rowArray_[0]->getIndices());
919	dubiousWeights = matrix_->dubiousWeights(this, rowArray_[0]->getIndices());
920	#endif
921	// If values pass then get list of candidates
922	int * candidateList = NULL;
923	int numberCandidates = 0;
924	#ifdef CLP_DEBUG
925	bool wasInValuesPass = (givenDuals != NULL);
926	#endif
927	int candidate = -1;
928	if (givenDuals) {
929	assert (ifValuesPass);
930	ifValuesPass = 1;
931	candidateList = new int[numberRows_];
932	// move reduced costs across
933	CoinMemcpyN(givenDuals, numberRows_ + numberColumns_, dj_);
934	int iRow;
935	for (iRow = 0; iRow < numberRows_; iRow++) {
936	int iPivot = pivotVariable_[iRow];
937	if (flagged(iPivot))
938	continue;
939	if (fabs(dj_[iPivot]) > dualTolerance_) {
940	// for now safer to ignore free ones
941	if (lower_[iPivot] > -1.0e50 \|\| upper_[iPivot] < 1.0e50)
942	if (pivoted(iPivot))
943	candidateList[numberCandidates++] = iRow;
944	} else {
945	clearPivoted(iPivot);
946	}
947	}
948	// and set first candidate
949	if (!numberCandidates) {
950	delete [] candidateList;
951	delete [] givenDuals;
952	givenDuals = NULL;
953	candidateList = NULL;
954	int iRow;
955	for (iRow = 0; iRow < numberRows_; iRow++) {
956	int iPivot = pivotVariable_[iRow];
957	clearPivoted(iPivot);
958	}
959	}
960	} else {
961	assert (!ifValuesPass);
962	}
963	#ifdef CHECK_ACCURACY
964	{
965	if (numberIterations_) {
966	int il = -1;
967	double largest = 1.0e-1;
968	int ilnb = -1;
969	double largestnb = 1.0e-8;
970	for (int i = 0; i < numberRows_ + numberColumns_; i++) {
971	double diff = fabs(solution_[i] - zzzzzz[i]);
972	if (diff > largest) {
973	largest = diff;
974	il = i;
975	}
976	if (getColumnStatus(i) != basic) {
977	if (diff > largestnb) {
978	largestnb = diff;
979	ilnb = i;
980	}
981	}
982	}
983	if (il >= 0 && ilnb < 0)
984	printf("largest diff of %g at %d, nonbasic %g at %d\n",
985	largest, il, largestnb, ilnb);
986	}
987	}
988	#endif
989	while (problemStatus_ == -1) {
990	//if (numberIterations_>=101624)
991	//resetFakeBounds(-1);
992	#ifdef CLP_DEBUG
993	if (givenDuals) {
994	double value5 = 0.0;
995	int i;
996	for (i = 0; i < numberRows_ + numberColumns_; i++) {
997	if (dj_[i] < -1.0e-6)
998	if (upper_[i] < 1.0e20)
999	value5 += dj_[i] * upper_[i];
1000	else
1001	printf("bad dj %g on %d with large upper status %d\n",
1002	dj_[i], i, status_[i] & 7);
1003	else if (dj_[i] > 1.0e-6)
1004	if (lower_[i] > -1.0e20)
1005	value5 += dj_[i] * lower_[i];
1006	else
1007	printf("bad dj %g on %d with large lower status %d\n",
1008	dj_[i], i, status_[i] & 7);
1009	}
1010	printf("Values objective Value %g\n", value5);
1011	}
1012	if ((handler_->logLevel() & 32) && wasInValuesPass) {
1013	double value5 = 0.0;
1014	int i;
1015	for (i = 0; i < numberRows_ + numberColumns_; i++) {
1016	if (dj_[i] < -1.0e-6)
1017	if (upper_[i] < 1.0e20)
1018	value5 += dj_[i] * upper_[i];
1019	else if (dj_[i] > 1.0e-6)
1020	if (lower_[i] > -1.0e20)
1021	value5 += dj_[i] * lower_[i];
1022	}
1023	printf("Values objective Value %g\n", value5);
1024	{
1025	int i;
1026	for (i = 0; i < numberRows_ + numberColumns_; i++) {
1027	int iSequence = i;
1028	double oldValue;
1029
1030	switch(getStatus(iSequence)) {
1031
1032	case basic:
1033	case ClpSimplex::isFixed:
1034	break;
1035	case isFree:
1036	case superBasic:
1037	abort();
1038	break;
1039	case atUpperBound:
1040	oldValue = dj_[iSequence];
1041	//assert (oldValue<=tolerance);
1042	assert (fabs(solution_[iSequence] - upper_[iSequence]) < 1.0e-7);
1043	break;
1044	case atLowerBound:
1045	oldValue = dj_[iSequence];
1046	//assert (oldValue>=-tolerance);
1047	assert (fabs(solution_[iSequence] - lower_[iSequence]) < 1.0e-7);
1048	break;
1049	}
1050	}
1051	}
1052	}
1053	#endif
1054	#ifdef CLP_DEBUG
1055	{
1056	int i;
1057	for (i = 0; i < 4; i++) {
1058	rowArray_[i]->checkClear();
1059	}
1060	for (i = 0; i < 2; i++) {
1061	columnArray_[i]->checkClear();
1062	}
1063	}
1064	#endif
1065	#if CLP_DEBUG>2
1066	// very expensive
1067	if (numberIterations_ > 3063 && numberIterations_ < 30700) {
1068	//handler_->setLogLevel(63);
1069	double saveValue = objectiveValue_;
1070	double * saveRow1 = new double[numberRows_];
1071	double * saveRow2 = new double[numberRows_];
1072	CoinMemcpyN(rowReducedCost_, numberRows_, saveRow1);
1073	CoinMemcpyN(rowActivityWork_, numberRows_, saveRow2);
1074	double * saveColumn1 = new double[numberColumns_];
1075	double * saveColumn2 = new double[numberColumns_];
1076	CoinMemcpyN(reducedCostWork_, numberColumns_, saveColumn1);
1077	CoinMemcpyN(columnActivityWork_, numberColumns_, saveColumn2);
1078	gutsOfSolution(NULL, NULL);
1079	printf("xxx %d old obj %g, recomputed %g, sum dual inf %g\n",
1080	numberIterations_,
1081	saveValue, objectiveValue_, sumDualInfeasibilities_);
1082	if (saveValue > objectiveValue_ + 1.0e-2)
1083	printf("bad\n");
1084	CoinMemcpyN(saveRow1, numberRows_, rowReducedCost_);
1085	CoinMemcpyN(saveRow2, numberRows_, rowActivityWork_);
1086	CoinMemcpyN(saveColumn1, numberColumns_, reducedCostWork_);
1087	CoinMemcpyN(saveColumn2, numberColumns_, columnActivityWork_);
1088	delete [] saveRow1;
1089	delete [] saveRow2;
1090	delete [] saveColumn1;
1091	delete [] saveColumn2;
1092	objectiveValue_ = saveValue;
1093	}
1094	#endif
1095	#if 0
1096	// if (factorization_->pivots()){
1097	{
1098	int iPivot;
1099	double * array = rowArray_[3]->denseVector();
1100	int i;
1101	for (iPivot = 0; iPivot < numberRows_; iPivot++) {
1102	int iSequence = pivotVariable_[iPivot];
1103	unpack(rowArray_[3], iSequence);
1104	factorization_->updateColumn(rowArray_[2], rowArray_[3]);
1105	assert (fabs(array[iPivot] - 1.0) < 1.0e-4);
1106	array[iPivot] = 0.0;
1107	for (i = 0; i < numberRows_; i++)
1108	assert (fabs(array[i]) < 1.0e-4);
1109	rowArray_[3]->clear();
1110	}
1111	}
1112	#endif
1113	#ifdef CLP_DEBUG
1114	{
1115	int iSequence, number = numberRows_ + numberColumns_;
1116	for (iSequence = 0; iSequence < number; iSequence++) {
1117	double lowerValue = lower_[iSequence];
1118	double upperValue = upper_[iSequence];
1119	double value = solution_[iSequence];
1120	if(getStatus(iSequence) != basic && getStatus(iSequence) != isFree) {
1121	assert(lowerValue > -1.0e20);
1122	assert(upperValue < 1.0e20);
1123	}
1124	switch(getStatus(iSequence)) {
1125
1126	case basic:
1127	break;
1128	case isFree:
1129	case superBasic:
1130	break;
1131	case atUpperBound:
1132	assert (fabs(value - upperValue) <= primalTolerance_) ;
1133	break;
1134	case atLowerBound:
1135	case ClpSimplex::isFixed:
1136	assert (fabs(value - lowerValue) <= primalTolerance_) ;
1137	break;
1138	}
1139	}
1140	}
1141	if(numberIterations_ == debugIteration) {
1142	printf("dodgy iteration coming up\n");
1143	}
1144	#endif
1145	#if 0
1146	printf("checking nz\n");
1147	for (int i = 0; i < 3; i++) {
1148	if (!rowArray_[i]->getNumElements())
1149	rowArray_[i]->checkClear();
1150	}
1151	#endif
1152	// choose row to go out
1153	// dualRow will go to virtual row pivot choice algorithm
1154	// make sure values pass off if it should be
1155	if (numberCandidates)
1156	candidate = candidateList[--numberCandidates];
1157	else
1158	candidate = -1;
1159	dualRow(candidate);
1160	if (pivotRow_ >= 0) {
1161	// we found a pivot row
1162	if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
1163	handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
1164	<< pivotRow_
1165	<< CoinMessageEol;
1166	}
1167	// check accuracy of weights
1168	dualRowPivot_->checkAccuracy();
1169	// Get good size for pivot
1170	// Allow first few iterations to take tiny
1171	double acceptablePivot = 1.0e-1 * acceptablePivot_;
1172	if (numberIterations_ > 100)
1173	acceptablePivot = acceptablePivot_;
1174	if (factorization_->pivots() > 10 \|\|
1175	(factorization_->pivots() && saveSumDual))
1176	acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
1177	else if (factorization_->pivots() > 5)
1178	acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
1179	else if (factorization_->pivots())
1180	acceptablePivot = acceptablePivot_; // relax
1181	// But factorizations complain if <1.0e-8
1182	//acceptablePivot=CoinMax(acceptablePivot,1.0e-8);
1183	double bestPossiblePivot = 1.0;
1184	// get sign for finding row of tableau
1185	if (candidate < 0) {
1186	// normal iteration
1187	// create as packed
1188	double direction = directionOut_;
1189	rowArray_[0]->createPacked(1, &pivotRow_, &direction);
1190	factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
1191	// Allow to do dualColumn0
1192	if (numberThreads_ < -1)
1193	spareIntArray_[0] = 1;
1194	spareDoubleArray_[0] = acceptablePivot;
1195	rowArray_[3]->clear();
1196	sequenceIn_ = -1;
1197	// put row of tableau in rowArray[0] and columnArray[0]
1198	assert (!rowArray_[1]->getNumElements());
1199	if (!scaledMatrix_) {
1200	if ((moreSpecialOptions_ & 8) != 0 && !rowScale_)
1201	spareIntArray_[0] = 1;
1202	matrix_->transposeTimes(this, -1.0,
1203	rowArray_[0], rowArray_[1], columnArray_[0]);
1204	} else {
1205	double * saveR = rowScale_;
1206	double * saveC = columnScale_;
1207	rowScale_ = NULL;
1208	columnScale_ = NULL;
1209	if ((moreSpecialOptions_ & 8) != 0)
1210	spareIntArray_[0] = 1;
1211	scaledMatrix_->transposeTimes(this, -1.0,
1212	rowArray_[0], rowArray_[1], columnArray_[0]);
1213	rowScale_ = saveR;
1214	columnScale_ = saveC;
1215	}
1216	#ifdef CLP_REPORT_PROGRESS
1217	memcpy(savePSol, solution_, (numberColumns_ + numberRows_)*sizeof(double));
1218	memcpy(saveDj, dj_, (numberColumns_ + numberRows_)*sizeof(double));
1219	memcpy(saveCost, cost_, (numberColumns_ + numberRows_)*sizeof(double));
1220	memcpy(saveStat, status_, (numberColumns_ + numberRows_)*sizeof(char));
1221	#endif
1222	// do ratio test for normal iteration
1223	bestPossiblePivot = dualColumn(rowArray_[0], columnArray_[0], rowArray_[3],
1224	columnArray_[1], acceptablePivot, dubiousWeights);
1225	} else {
1226	// Make sure direction plausible
1227	CoinAssert (upperOut_ < 1.0e50 \|\| lowerOut_ > -1.0e50);
1228	// If in integer cleanup do direction using duals
1229	// may be wrong way round
1230	if(ifValuesPass == 2) {
1231	if (dual_[pivotRow_] > 0.0) {
1232	// this will give a -1 in pivot row (as slacks are -1.0)
1233	directionOut_ = 1;
1234	} else {
1235	directionOut_ = -1;
1236	}
1237	}
1238	if (directionOut_ < 0 && fabs(valueOut_ - upperOut_) > dualBound_ + primalTolerance_) {
1239	if (fabs(valueOut_ - upperOut_) > fabs(valueOut_ - lowerOut_))
1240	directionOut_ = 1;
1241	} else if (directionOut_ > 0 && fabs(valueOut_ - lowerOut_) > dualBound_ + primalTolerance_) {
1242	if (fabs(valueOut_ - upperOut_) < fabs(valueOut_ - lowerOut_))
1243	directionOut_ = -1;
1244	}
1245	double direction = directionOut_;
1246	rowArray_[0]->createPacked(1, &pivotRow_, &direction);
1247	factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
1248	// put row of tableau in rowArray[0] and columnArray[0]
1249	if (!scaledMatrix_) {
1250	matrix_->transposeTimes(this, -1.0,
1251	rowArray_[0], rowArray_[3], columnArray_[0]);
1252	} else {
1253	double * saveR = rowScale_;
1254	double * saveC = columnScale_;
1255	rowScale_ = NULL;
1256	columnScale_ = NULL;
1257	scaledMatrix_->transposeTimes(this, -1.0,
1258	rowArray_[0], rowArray_[3], columnArray_[0]);
1259	rowScale_ = saveR;
1260	columnScale_ = saveC;
1261	}
1262	acceptablePivot *= 10.0;
1263	// do ratio test
1264	if (ifValuesPass == 1) {
1265	checkPossibleValuesMove(rowArray_[0], columnArray_[0],
1266	acceptablePivot);
1267	} else {
1268	checkPossibleCleanup(rowArray_[0], columnArray_[0],
1269	acceptablePivot);
1270	if (sequenceIn_ < 0) {
1271	rowArray_[0]->clear();
1272	columnArray_[0]->clear();
1273	continue; // can't do anything
1274	}
1275	}
1276
1277	// recompute true dualOut_
1278	if (directionOut_ < 0) {
1279	dualOut_ = valueOut_ - upperOut_;
1280	} else {
1281	dualOut_ = lowerOut_ - valueOut_;
1282	}
1283	// check what happened if was values pass
1284	// may want to move part way i.e. movement
1285	bool normalIteration = (sequenceIn_ != sequenceOut_);
1286
1287	clearPivoted(sequenceOut_); // make sure won't be done again
1288	// see if end of values pass
1289	if (!numberCandidates) {
1290	int iRow;
1291	delete [] candidateList;
1292	delete [] givenDuals;
1293	candidate = -2; // -2 signals end
1294	givenDuals = NULL;
1295	candidateList = NULL;
1296	ifValuesPass = 1;
1297	for (iRow = 0; iRow < numberRows_; iRow++) {
1298	int iPivot = pivotVariable_[iRow];
1299	//assert (fabs(dj_[iPivot]),1.0e-5);
1300	clearPivoted(iPivot);
1301	}
1302	}
1303	if (!normalIteration) {
1304	//rowArray_[0]->cleanAndPackSafe(1.0e-60);
1305	//columnArray_[0]->cleanAndPackSafe(1.0e-60);
1306	updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
1307	if (candidate == -2)
1308	problemStatus_ = -2;
1309	continue; // skip rest of iteration
1310	} else {
1311	// recompute dualOut_
1312	if (directionOut_ < 0) {
1313	dualOut_ = valueOut_ - upperOut_;
1314	} else {
1315	dualOut_ = lowerOut_ - valueOut_;
1316	}
1317	}
1318	}
1319	if (sequenceIn_ >= 0) {
1320	// normal iteration
1321	// update the incoming column
1322	double btranAlpha = -alpha_ * directionOut_; // for check
1323	unpackPacked(rowArray_[1]);
1324	// moved into updateWeights - factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
1325	// and update dual weights (can do in parallel - with extra array)
1326	alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
1327	rowArray_[2],
1328	rowArray_[3],
1329	rowArray_[1]);
1330	// see if update stable
1331	#ifdef CLP_DEBUG
1332	if ((handler_->logLevel() & 32))
1333	printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
1334	#endif
1335	double checkValue = 1.0e-7;
1336	// if can't trust much and long way from optimal then relax
1337	if (largestPrimalError_ > 10.0)
1338	checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
1339	if (fabs(btranAlpha) < 1.0e-12 \|\| fabs(alpha_) < 1.0e-12 \|\|
1340	fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
1341	handler_->message(CLP_DUAL_CHECK, messages_)
1342	<< btranAlpha
1343	<< alpha_
1344	<< CoinMessageEol;
1345	if (factorization_->pivots()) {
1346	dualRowPivot_->unrollWeights();
1347	problemStatus_ = -2; // factorize now
1348	rowArray_[0]->clear();
1349	rowArray_[1]->clear();
1350	columnArray_[0]->clear();
1351	returnCode = -2;
1352	break;
1353	} else {
1354	// take on more relaxed criterion
1355	double test;
1356	if (fabs(btranAlpha) < 1.0e-8 \|\| fabs(alpha_) < 1.0e-8)
1357	test = 1.0e-1 * fabs(alpha_);
1358	else
1359	test = 1.0e-4 * (1.0 + fabs(alpha_));
1360	if (fabs(btranAlpha) < 1.0e-12 \|\| fabs(alpha_) < 1.0e-12 \|\|
1361	fabs(btranAlpha - alpha_) > test) {
1362	dualRowPivot_->unrollWeights();
1363	// need to reject something
1364	char x = isColumn(sequenceOut_) ? 'C' : 'R';
1365	handler_->message(CLP_SIMPLEX_FLAG, messages_)
1366	<< x << sequenceWithin(sequenceOut_)
1367	<< CoinMessageEol;
1368	#ifdef COIN_DEVELOP
1369	printf("flag a %g %g\n", btranAlpha, alpha_);
1370	#endif
1371	//#define FEB_TRY
1372	#if 1 //def FEB_TRY
1373	// Make safer?
1374	factorization_->saferTolerances (-0.99, -1.03);
1375	#endif
1376	setFlagged(sequenceOut_);
1377	progress_.clearBadTimes();
1378	lastBadIteration_ = numberIterations_; // say be more cautious
1379	rowArray_[0]->clear();
1380	rowArray_[1]->clear();
1381	columnArray_[0]->clear();
1382	if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
1383	//printf("I think should declare infeasible\n");
1384	problemStatus_ = 1;
1385	returnCode = 1;
1386	break;
1387	}
1388	continue;
1389	}
1390	}
1391	}
1392	// update duals BEFORE replaceColumn so can do updateColumn
1393	double objectiveChange = 0.0;
1394	// do duals first as variables may flip bounds
1395	// rowArray_[0] and columnArray_[0] may have flips
1396	// so use rowArray_[3] for work array from here on
1397	int nswapped = 0;
1398	//rowArray_[0]->cleanAndPackSafe(1.0e-60);
1399	//columnArray_[0]->cleanAndPackSafe(1.0e-60);
1400	if (candidate == -1) {
1401	// make sure incoming doesn't count
1402	Status saveStatus = getStatus(sequenceIn_);
1403	setStatus(sequenceIn_, basic);
1404	nswapped = updateDualsInDual(rowArray_[0], columnArray_[0],
1405	rowArray_[2], theta_,
1406	objectiveChange, false);
1407	setStatus(sequenceIn_, saveStatus);
1408	} else {
1409	updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
1410	}
1411	double oldDualOut = dualOut_;
1412	// which will change basic solution
1413	if (nswapped) {
1414	if (rowArray_[2]->getNumElements()) {
1415	factorization_->updateColumn(rowArray_[3], rowArray_[2]);
1416	dualRowPivot_->updatePrimalSolution(rowArray_[2],
1417	1.0, objectiveChange);
1418	}
1419	// recompute dualOut_
1420	valueOut_ = solution_[sequenceOut_];
1421	if (directionOut_ < 0) {
1422	dualOut_ = valueOut_ - upperOut_;
1423	} else {
1424	dualOut_ = lowerOut_ - valueOut_;
1425	}
1426	#if 0
1427	if (dualOut_ < 0.0) {
1428	#ifdef CLP_DEBUG
1429	if (handler_->logLevel() & 32) {
1430	printf(" dualOut_ %g %g save %g\n", dualOut_, averagePrimalInfeasibility, saveDualOut);
1431	printf("values %g %g %g %g %g %g %g\n", lowerOut_, valueOut_, upperOut_,
1432	objectiveChange,);
1433	}
1434	#endif
1435	if (upperOut_ == lowerOut_)
1436	dualOut_ = 0.0;
1437	}
1438	if(dualOut_ < -CoinMax(1.0e-12 * averagePrimalInfeasibility, 1.0e-8)
1439	&& factorization_->pivots() > 100 &&
1440	getStatus(sequenceIn_) != isFree) {
1441	// going backwards - factorize
1442	dualRowPivot_->unrollWeights();
1443	problemStatus_ = -2; // factorize now
1444	returnCode = -2;
1445	break;
1446	}
1447	#endif
1448	}
1449	// amount primal will move
1450	double movement = -dualOut_ * directionOut_ / alpha_;
1451	double movementOld = oldDualOut * directionOut_ / alpha_;
1452	// so objective should increase by fabs(dj)*movement
1453	// but we already have objective change - so check will be good
1454	if (objectiveChange + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(objectiveValue_))) {
1455	#ifdef CLP_DEBUG
1456	if (handler_->logLevel() & 32)
1457	printf("movement %g, swap change %g, rest %g * %g\n",
1458	objectiveChange + fabs(movement * dualIn_),
1459	objectiveChange, movement, dualIn_);
1460	#endif
1461	if(factorization_->pivots()) {
1462	// going backwards - factorize
1463	dualRowPivot_->unrollWeights();
1464	problemStatus_ = -2; // factorize now
1465	returnCode = -2;
1466	break;
1467	}
1468	}
1469	// if stable replace in basis
1470	int updateStatus = factorization_->replaceColumn(this,
1471	rowArray_[2],
1472	rowArray_[1],
1473	pivotRow_,
1474	alpha_,
1475	(moreSpecialOptions_ & 16) != 0,
1476	acceptablePivot);
1477	// If looks like bad pivot - refactorize
1478	if (fabs(dualOut_) > 1.0e50)
1479	updateStatus = 2;
1480	// if no pivots, bad update but reasonable alpha - take and invert
1481	if (updateStatus == 2 &&
1482	!factorization_->pivots() && fabs(alpha_) > 1.0e-5)
1483	updateStatus = 4;
1484	if (updateStatus == 1 \|\| updateStatus == 4) {
1485	// slight error
1486	if (factorization_->pivots() > 5 \|\| updateStatus == 4) {
1487	problemStatus_ = -2; // factorize now
1488	returnCode = -3;
1489	}
1490	} else if (updateStatus == 2) {
1491	// major error
1492	dualRowPivot_->unrollWeights();
1493	// later we may need to unwind more e.g. fake bounds
1494	if (factorization_->pivots() &&
1495	((moreSpecialOptions_ & 16) == 0 \|\| factorization_->pivots() > 4)) {
1496	problemStatus_ = -2; // factorize now
1497	returnCode = -2;
1498	moreSpecialOptions_ \|= 16;
1499	break;
1500	} else {
1501	// need to reject something
1502	char x = isColumn(sequenceOut_) ? 'C' : 'R';
1503	handler_->message(CLP_SIMPLEX_FLAG, messages_)
1504	<< x << sequenceWithin(sequenceOut_)
1505	<< CoinMessageEol;
1506	#ifdef COIN_DEVELOP
1507	printf("flag b %g\n", alpha_);
1508	#endif
1509	setFlagged(sequenceOut_);
1510	progress_.clearBadTimes();
1511	lastBadIteration_ = numberIterations_; // say be more cautious
1512	rowArray_[0]->clear();
1513	rowArray_[1]->clear();
1514	columnArray_[0]->clear();
1515	// make sure dual feasible
1516	// look at all rows and columns
1517	double objectiveChange = 0.0;
1518	updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
1519	0.0, objectiveChange, true);
1520	rowArray_[1]->clear();
1521	columnArray_[0]->clear();
1522	continue;
1523	}
1524	} else if (updateStatus == 3) {
1525	// out of memory
1526	// increase space if not many iterations
1527	if (factorization_->pivots() <
1528	0.5 * factorization_->maximumPivots() &&
1529	factorization_->pivots() < 200)
1530	factorization_->areaFactor(
1531	factorization_->areaFactor() * 1.1);
1532	problemStatus_ = -2; // factorize now
1533	} else if (updateStatus == 5) {
1534	problemStatus_ = -2; // factorize now
1535	}
1536	// update primal solution
1537	if (theta_ < 0.0 && candidate == -1) {
1538	#ifdef CLP_DEBUG
1539	if (handler_->logLevel() & 32)
1540	printf("negative theta %g\n", theta_);
1541	#endif
1542	theta_ = 0.0;
1543	}
1544	// do actual flips
1545	flipBounds(rowArray_[0], columnArray_[0]);
1546	//rowArray_[1]->expand();
1547	dualRowPivot_->updatePrimalSolution(rowArray_[1],
1548	movement,
1549	objectiveChange);
1550	#ifdef CLP_DEBUG
1551	double oldobj = objectiveValue_;
1552	#endif
1553	// modify dualout
1554	dualOut_ /= alpha_;
1555	dualOut_ *= -directionOut_;
1556	//setStatus(sequenceIn_,basic);
1557	dj_[sequenceIn_] = 0.0;
1558	double oldValue = valueIn_;
1559	if (directionIn_ == -1) {
1560	// as if from upper bound
1561	valueIn_ = upperIn_ + dualOut_;
1562	} else {
1563	// as if from lower bound
1564	valueIn_ = lowerIn_ + dualOut_;
1565	}
1566	objectiveChange += cost_[sequenceIn_] * (valueIn_ - oldValue);
1567	// outgoing
1568	// set dj to zero unless values pass
1569	if (directionOut_ > 0) {
1570	valueOut_ = lowerOut_;
1571	if (candidate == -1)
1572	dj_[sequenceOut_] = theta_;
1573	} else {
1574	valueOut_ = upperOut_;
1575	if (candidate == -1)
1576	dj_[sequenceOut_] = -theta_;
1577	}
1578	solution_[sequenceOut_] = valueOut_;
1579	int whatNext = housekeeping(objectiveChange);
1580	#ifdef CLP_REPORT_PROGRESS
1581	if (ixxxxxx > ixxyyyy - 5) {
1582	handler_->setLogLevel(63);
1583	int nTotal = numberColumns_ + numberRows_;
1584	double oldObj = 0.0;
1585	double newObj = 0.0;
1586	for (int i = 0; i < nTotal; i++) {
1587	if (savePSol[i])
1588	oldObj += savePSol[i] * saveCost[i];
1589	if (solution_[i])
1590	newObj += solution_[i] * cost_[i];
1591	bool printIt = false;
1592	if (cost_[i] != saveCost[i])
1593	printIt = true;
1594	if (status_[i] != saveStat[i])
1595	printIt = true;
1596	if (printIt)
1597	printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
1598	i, saveStat[i], saveCost[i], savePSol[i],
1599	status_[i], cost_[i], solution_[i]);
1600	// difference
1601	savePSol[i] = solution_[i] - savePSol[i];
1602	}
1603	printf("pivots %d, old obj %g new %g\n",
1604	factorization_->pivots(),
1605	oldObj, newObj);
1606	memset(saveDj, 0, numberRows_ * sizeof(double));
1607	times(1.0, savePSol, saveDj);
1608	double largest = 1.0e-6;
1609	int k = -1;
1610	for (int i = 0; i < numberRows_; i++) {
1611	saveDj[i] -= savePSol[i+numberColumns_];
1612	if (fabs(saveDj[i]) > largest) {
1613	largest = fabs(saveDj[i]);
1614	k = i;
1615	}
1616	}
1617	if (k >= 0)
1618	printf("Not null %d %g\n", k, largest);
1619	}
1620	#endif
1621	#ifdef VUB
1622	{
1623	if ((sequenceIn_ < numberColumns_ && vub[sequenceIn_] >= 0) \|\| toVub[sequenceIn_] >= 0 \|\|
1624	(sequenceOut_ < numberColumns_ && vub[sequenceOut_] >= 0) \|\| toVub[sequenceOut_] >= 0) {
1625	int inSequence = sequenceIn_;
1626	int inVub = -1;
1627	if (sequenceIn_ < numberColumns_)
1628	inVub = vub[sequenceIn_];
1629	int inBack = toVub[inSequence];
1630	int inSlack = -1;
1631	if (inSequence >= numberColumns_ && inBack >= 0) {
1632	inSlack = inSequence - numberColumns_;
1633	inSequence = inBack;
1634	inBack = toVub[inSequence];
1635	}
1636	if (inVub >= 0)
1637	printf("Vub %d in ", inSequence);
1638	if (inBack >= 0 && inSlack < 0)
1639	printf("%d (descendent of %d) in ", inSequence, inBack);
1640	if (inSlack >= 0)
1641	printf("slack for row %d -> %d (descendent of %d) in ", inSlack, inSequence, inBack);
1642	int outSequence = sequenceOut_;
1643	int outVub = -1;
1644	if (sequenceOut_ < numberColumns_)
1645	outVub = vub[sequenceOut_];
1646	int outBack = toVub[outSequence];
1647	int outSlack = -1;
1648	if (outSequence >= numberColumns_ && outBack >= 0) {
1649	outSlack = outSequence - numberColumns_;
1650	outSequence = outBack;
1651	outBack = toVub[outSequence];
1652	}
1653	if (outVub >= 0)
1654	printf("Vub %d out ", outSequence);
1655	if (outBack >= 0 && outSlack < 0)
1656	printf("%d (descendent of %d) out ", outSequence, outBack);
1657	if (outSlack >= 0)
1658	printf("slack for row %d -> %d (descendent of %d) out ", outSlack, outSequence, outBack);
1659	printf("\n");
1660	}
1661	}
1662	#endif
1663	#if 0
1664	if (numberIterations_ > 206033)
1665	handler_->setLogLevel(63);
1666	if (numberIterations_ > 210567)
1667	exit(77);
1668	#endif
1669	if (!givenDuals && ifValuesPass && ifValuesPass != 2) {
1670	handler_->message(CLP_END_VALUES_PASS, messages_)
1671	<< numberIterations_;
1672	whatNext = 1;
1673	}
1674	#ifdef CHECK_ACCURACY
1675	if (whatNext) {
1676	CoinMemcpyN(solution_, (numberRows_ + numberColumns_), zzzzzz);
1677	}
1678	#endif
1679	//if (numberIterations_==1890)
1680	//whatNext=1;
1681	//if (numberIterations_>2000)
1682	//exit(77);
1683	// and set bounds correctly
1684	originalBound(sequenceIn_);
1685	changeBound(sequenceOut_);
1686	#ifdef CLP_DEBUG
1687	if (objectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16))
1688	printf("obj backwards %g %g\n", objectiveValue_, oldobj);
1689	#endif
1690	#if 0
1691	{
1692	for (int i = 0; i < numberRows_ + numberColumns_; i++) {
1693	FakeBound bound = getFakeBound(i);
1694	if (bound == ClpSimplexDual::upperFake) {
1695	assert (upper_[i] < 1.0e20);
1696	} else if (bound == ClpSimplexDual::lowerFake) {
1697	assert (lower_[i] > -1.0e20);
1698	} else if (bound == ClpSimplexDual::bothFake) {
1699	assert (upper_[i] < 1.0e20);
1700	assert (lower_[i] > -1.0e20);
1701	}
1702	}
1703	}
1704	#endif
1705	if (whatNext == 1 \|\| candidate == -2) {
1706	problemStatus_ = -2; // refactorize
1707	} else if (whatNext == 2) {
1708	// maximum iterations or equivalent
1709	problemStatus_ = 3;
1710	returnCode = 3;
1711	break;
1712	}
1713	// Check event
1714	{
1715	int status = eventHandler_->event(ClpEventHandler::endOfIteration);
1716	if (status >= 0) {
1717	problemStatus_ = 5;
1718	secondaryStatus_ = ClpEventHandler::endOfIteration;
1719	returnCode = 4;
1720	break;
1721	}
1722	}
1723	} else {
1724	#ifdef CLP_INVESTIGATE_SERIAL
1725	z_thinks = 1;
1726	#endif
1727	// no incoming column is valid
1728	pivotRow_ = -1;
1729	#ifdef CLP_DEBUG
1730	if (handler_->logLevel() & 32)
1731	printf("** no column pivot\n");
1732	#endif
1733	if (factorization_->pivots() < 2 && acceptablePivot_ <= 1.0e-8) {
1734	//&&goodAccuracy()) {
1735	// If not in branch and bound etc save ray
1736	delete [] ray_;
1737	if ((specialOptions_&(1024 \| 4096)) == 0 \|\| (specialOptions_ & 32) != 0) {
1738	// create ray anyway
1739	ray_ = new double [ numberRows_];
1740	rowArray_[0]->expand(); // in case packed
1741	CoinMemcpyN(rowArray_[0]->denseVector(), numberRows_, ray_);
1742	} else {
1743	ray_ = NULL;
1744	}
1745	// If we have just factorized and infeasibility reasonable say infeas
1746	double dualTest = ((specialOptions_ & 4096) != 0) ? 1.0e8 : 1.0e13;
1747	if (((specialOptions_ & 4096) != 0 \|\| bestPossiblePivot < 1.0e-11) && dualBound_ > dualTest) {
1748	double testValue = 1.0e-4;
1749	if (!factorization_->pivots() && numberPrimalInfeasibilities_ == 1)
1750	testValue = 1.0e-6;
1751	if (valueOut_ > upperOut_ + testValue \|\| valueOut_ < lowerOut_ - testValue
1752	\|\| (specialOptions_ & 64) == 0) {
1753	// say infeasible
1754	problemStatus_ = 1;
1755	// unless primal feasible!!!!
1756	//printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
1757	// numberDualInfeasibilities_,sumDualInfeasibilities_);
1758	//#define TEST_CLP_NODE
1759	#ifndef TEST_CLP_NODE
1760	// Should be correct - but ...
1761	int numberFake = numberAtFakeBound();
1762	double sumPrimal = (!numberFake) ? 2.0e5 : sumPrimalInfeasibilities_;
1763	if (sumPrimalInfeasibilities_ < 1.0e-3 \|\| sumDualInfeasibilities_ > 1.0e-5 \|\|
1764	(sumPrimal < 1.0e5 && (specialOptions_ & 1024) != 0 && factorization_->pivots())) {
1765	if (sumPrimal > 50.0 && factorization_->pivots() > 2) {
1766	problemStatus_ = -4;
1767	#ifdef COIN_DEVELOP
1768	printf("status to -4 at %d - primalinf %g pivots %d\n",
1769	__LINE__, sumPrimalInfeasibilities_,
1770	factorization_->pivots());
1771	#endif
1772	} else {
1773	problemStatus_ = 10;
1774	#if COIN_DEVELOP>1
1775	printf("returning at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d - options (1024-16384) %d %d %d %d %d\n",
1776	__LINE__, numberPrimalInfeasibilities_,
1777	sumPrimalInfeasibilities_,
1778	numberDualInfeasibilities_, sumDualInfeasibilities_,
1779	numberFake_, dualBound_, factorization_->pivots(),
1780	(specialOptions_ & 1024) != 0 ? 1 : 0,
1781	(specialOptions_ & 2048) != 0 ? 1 : 0,
1782	(specialOptions_ & 4096) != 0 ? 1 : 0,
1783	(specialOptions_ & 8192) != 0 ? 1 : 0,
1784	(specialOptions_ & 16384) != 0 ? 1 : 0
1785	);
1786	#endif
1787	// Get rid of objective
1788	if ((specialOptions_ & 16384) == 0)
1789	objective_ = new ClpLinearObjective(NULL, numberColumns_);
1790	}
1791	}
1792	#else
1793	if (sumPrimalInfeasibilities_ < 1.0e-3 \|\| sumDualInfeasibilities_ > 1.0e-6) {
1794	#ifdef COIN_DEVELOP
1795	printf("at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d\n",
1796	__LINE__, numberPrimalInfeasibilities_,
1797	sumPrimalInfeasibilities_,
1798	numberDualInfeasibilities_, sumDualInfeasibilities_,
1799	numberFake_, dualBound_, factorization_->pivots());
1800	#endif
1801	if ((specialOptions_ & 1024) != 0 && factorization_->pivots()) {
1802	problemStatus_ = 10;
1803	#if COIN_DEVELOP>1
1804	printf("returning at %d\n", __LINE__);
1805	#endif
1806	// Get rid of objective
1807	if ((specialOptions_ & 16384) == 0)
1808	objective_ = new ClpLinearObjective(NULL, numberColumns_);
1809	}
1810	}
1811	#endif
1812	rowArray_[0]->clear();
1813	columnArray_[0]->clear();
1814	returnCode = 1;
1815	break;
1816	}
1817	}
1818	// If special option set - put off as long as possible
1819	if ((specialOptions_ & 64) == 0 \|\| (moreSpecialOptions_ & 64) != 0) {
1820	if (factorization_->pivots() == 0)
1821	problemStatus_ = -4; //say looks infeasible
1822	} else {
1823	// flag
1824	char x = isColumn(sequenceOut_) ? 'C' : 'R';
1825	handler_->message(CLP_SIMPLEX_FLAG, messages_)
1826	<< x << sequenceWithin(sequenceOut_)
1827	<< CoinMessageEol;
1828	#ifdef COIN_DEVELOP
1829	printf("flag c\n");
1830	#endif
1831	setFlagged(sequenceOut_);
1832	if (!factorization_->pivots()) {
1833	rowArray_[0]->clear();
1834	columnArray_[0]->clear();
1835	continue;
1836	}
1837	}
1838	}
1839	if (factorization_->pivots() < 5 && acceptablePivot_ > 1.0e-8)
1840	acceptablePivot_ = 1.0e-8;
1841	rowArray_[0]->clear();
1842	columnArray_[0]->clear();
1843	returnCode = 1;
1844	break;
1845	}
1846	} else {
1847	#ifdef CLP_INVESTIGATE_SERIAL
1848	z_thinks = 0;
1849	#endif
1850	// no pivot row
1851	#ifdef CLP_DEBUG
1852	if (handler_->logLevel() & 32)
1853	printf("** no row pivot\n");
1854	#endif
1855	// If in branch and bound try and get rid of fixed variables
1856	if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
1857	assert (!candidateList);
1858	candidateList = new int[numberRows_];
1859	int iRow;
1860	for (iRow = 0; iRow < numberRows_; iRow++) {
1861	int iPivot = pivotVariable_[iRow];
1862	if (flagged(iPivot) \|\| !pivoted(iPivot))
1863	continue;
1864	assert (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]);
1865	candidateList[numberCandidates++] = iRow;
1866	}
1867	// and set first candidate
1868	if (!numberCandidates) {
1869	delete [] candidateList;
1870	candidateList = NULL;
1871	int iRow;
1872	for (iRow = 0; iRow < numberRows_; iRow++) {
1873	int iPivot = pivotVariable_[iRow];
1874	clearPivoted(iPivot);
1875	}
1876	} else {
1877	ifValuesPass = 2;
1878	continue;
1879	}
1880	}
1881	int numberPivots = factorization_->pivots();
1882	bool specialCase;
1883	int useNumberFake;
1884	returnCode = 0;
1885	if (numberPivots <= CoinMax(dontFactorizePivots_, 20) &&
1886	(specialOptions_ & 2048) != 0 && (true \|\| !numberChanged_ \|\| perturbation_ == 101)
1887	&& dualBound_ >= 1.0e8) {
1888	specialCase = true;
1889	// as dual bound high - should be okay
1890	useNumberFake = 0;
1891	} else {
1892	specialCase = false;
1893	useNumberFake = numberFake_;
1894	}
1895	if (!numberPivots \|\| specialCase) {
1896	// may have crept through - so may be optimal
1897	// check any flagged variables
1898	int iRow;
1899	for (iRow = 0; iRow < numberRows_; iRow++) {
1900	int iPivot = pivotVariable_[iRow];
1901	if (flagged(iPivot))
1902	break;
1903	}
1904	if (iRow < numberRows_ && numberPivots) {
1905	// try factorization
1906	returnCode = -2;
1907	}
1908
1909	if (useNumberFake \|\| numberDualInfeasibilities_) {
1910	// may be dual infeasible
1911	if ((specialOptions_ & 1024) == 0)
1912	problemStatus_ = -5;
1913	else if (!useNumberFake && numberPrimalInfeasibilities_
1914	&& !numberPivots)
1915	problemStatus_ = 1;
1916	} else {
1917	if (iRow < numberRows_) {
1918	#ifdef COIN_DEVELOP
1919	std::cout << "Flagged variables at end - infeasible?" << std::endl;
1920	printf("Probably infeasible - pivot was %g\n", alpha_);
1921	#endif
1922	//if (fabs(alpha_)<1.0e-4) {
1923	//problemStatus_=1;
1924	//} else {
1925	#ifdef CLP_DEBUG
1926	abort();
1927	#endif
1928	//}
1929	problemStatus_ = -5;
1930	} else {
1931	problemStatus_ = 0;
1932	#ifndef CLP_CHECK_NUMBER_PIVOTS
1933	#define CLP_CHECK_NUMBER_PIVOTS 10
1934	#endif
1935	#if CLP_CHECK_NUMBER_PIVOTS < 20
1936	if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
1937	#ifndef NDEBUG_CLP
1938	int nTotal = numberRows_ + numberColumns_;
1939	double * comp = CoinCopyOfArray(solution_, nTotal);
1940	#endif
1941	computePrimals(rowActivityWork_, columnActivityWork_);
1942	#ifndef NDEBUG_CLP
1943	double largest = 1.0e-5;
1944	int bad = -1;
1945	for (int i = 0; i < nTotal; i++) {
1946	double value = solution_[i];
1947	double larger = CoinMax(fabs(value), fabs(comp[i]));
1948	double tol = 1.0e-5 + 1.0e-5 * larger;
1949	double diff = fabs(value - comp[i]);
1950	if (diff - tol > largest) {
1951	bad = i;
1952	largest = diff - tol;
1953	}
1954	}
1955	if (bad >= 0)
1956	COIN_DETAIL_PRINT(printf("bad %d old %g new %g\n", bad, comp[bad], solution_[bad]));
1957	#endif
1958	checkPrimalSolution(rowActivityWork_, columnActivityWork_);
1959	if (numberPrimalInfeasibilities_) {
1960	#ifdef CLP_INVESTIGATE
1961	printf("XXX Infeas ? %d inf summing to %g\n", numberPrimalInfeasibilities_,
1962	sumPrimalInfeasibilities_);
1963	#endif
1964	problemStatus_ = -1;
1965	returnCode = -2;
1966	}
1967	#ifndef NDEBUG_CLP
1968	memcpy(solution_, comp, nTotal * sizeof(double));
1969	delete [] comp;
1970	#endif
1971	}
1972	#endif
1973	if (!problemStatus_) {
1974	// make it look OK
1975	numberPrimalInfeasibilities_ = 0;
1976	sumPrimalInfeasibilities_ = 0.0;
1977	numberDualInfeasibilities_ = 0;
1978	sumDualInfeasibilities_ = 0.0;
1979	// May be perturbed
1980	if (perturbation_ == 101 \|\| numberChanged_) {
1981	numberChanged_ = 0; // Number of variables with changed costs
1982	perturbation_ = 102; // stop any perturbations
1983	//double changeCost;
1984	//changeBounds(1,NULL,changeCost);
1985	createRim4(false);
1986	// make sure duals are current
1987	computeDuals(givenDuals);
1988	checkDualSolution();
1989	progress_.modifyObjective(-COIN_DBL_MAX);
1990	if (numberDualInfeasibilities_) {
1991	problemStatus_ = 10; // was -3;
1992	} else {
1993	computeObjectiveValue(true);
1994	}
1995	} else if (numberPivots) {
1996	computeObjectiveValue(true);
1997	}
1998	if (numberPivots < -1000) {
1999	// objective may be wrong
2000	objectiveValue_ = innerProduct(cost_, numberColumns_ + numberRows_, solution_);
2001	objectiveValue_ += objective_->nonlinearOffset();
2002	objectiveValue_ /= (objectiveScale_ * rhsScale_);
2003	if ((specialOptions_ & 16384) == 0) {
2004	// and dual_ may be wrong (i.e. for fixed or basic)
2005	CoinIndexedVector * arrayVector = rowArray_[1];
2006	arrayVector->clear();
2007	int iRow;
2008	double * array = arrayVector->denseVector();
2009	/* Use dual_ instead of array
2010	Even though dual_ is only numberRows_ long this is
2011	okay as gets permuted to longer rowArray_[2]
2012	*/
2013	arrayVector->setDenseVector(dual_);
2014	int * index = arrayVector->getIndices();
2015	int number = 0;
2016	for (iRow = 0; iRow < numberRows_; iRow++) {
2017	int iPivot = pivotVariable_[iRow];
2018	double value = cost_[iPivot];
2019	dual_[iRow] = value;
2020	if (value) {
2021	index[number++] = iRow;
2022	}
2023	}
2024	arrayVector->setNumElements(number);
2025	// Extended duals before "updateTranspose"
2026	matrix_->dualExpanded(this, arrayVector, NULL, 0);
2027	// Btran basic costs
2028	rowArray_[2]->clear();
2029	factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
2030	// and return vector
2031	arrayVector->setDenseVector(array);
2032	}
2033	}
2034	sumPrimalInfeasibilities_ = 0.0;
2035	}
2036	if ((specialOptions_&(1024 + 16384)) != 0 && !problemStatus_) {
2037	CoinIndexedVector * arrayVector = rowArray_[1];
2038	arrayVector->clear();
2039	double * rhs = arrayVector->denseVector();
2040	times(1.0, solution_, rhs);
2041	#ifdef CHECK_ACCURACY
2042	bool bad = false;
2043	#endif
2044	bool bad2 = false;
2045	int i;
2046	for ( i = 0; i < numberRows_; i++) {
2047	if (rhs[i] < rowLowerWork_[i] - primalTolerance_ \|\|
2048	rhs[i] > rowUpperWork_[i] + primalTolerance_) {
2049	bad2 = true;
2050	#ifdef CHECK_ACCURACY
2051	printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
2052	rowLowerWork_[i], rowUpperWork_[i],
2053	rhs[i], rowActivityWork_[i]);
2054	#endif
2055	} else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
2056	#ifdef CHECK_ACCURACY
2057	bad = true;
2058	printf("row %d correct %g bad %g\n", i, rhs[i], rowActivityWork_[i]);
2059	#endif
2060	}
2061	rhs[i] = 0.0;
2062	}
2063	for ( i = 0; i < numberColumns_; i++) {
2064	if (solution_[i] < columnLowerWork_[i] - primalTolerance_ \|\|
2065	solution_[i] > columnUpperWork_[i] + primalTolerance_) {
2066	bad2 = true;
2067	#ifdef CHECK_ACCURACY
2068	printf("column %d out of bounds %g, %g correct %g bad %g\n", i,
2069	columnLowerWork_[i], columnUpperWork_[i],
2070	solution_[i], columnActivityWork_[i]);
2071	#endif
2072	}
2073	}
2074	if (bad2) {
2075	problemStatus_ = -3;
2076	returnCode = -2;
2077	// Force to re-factorize early next time
2078	int numberPivots = factorization_->pivots();
2079	forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2080	}
2081	}
2082	}
2083	}
2084	} else {
2085	problemStatus_ = -3;
2086	returnCode = -2;
2087	// Force to re-factorize early next time
2088	int numberPivots = factorization_->pivots();
2089	forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2090	}
2091	break;
2092	}
2093	}
2094	if (givenDuals) {
2095	CoinMemcpyN(dj_, numberRows_ + numberColumns_, givenDuals);
2096	// get rid of any values pass array
2097	delete [] candidateList;
2098	}
2099	delete [] dubiousWeights;
2100	#ifdef CLP_REPORT_PROGRESS
2101	if (ixxxxxx > ixxyyyy - 5) {
2102	int nTotal = numberColumns_ + numberRows_;
2103	double oldObj = 0.0;
2104	double newObj = 0.0;
2105	for (int i = 0; i < nTotal; i++) {
2106	if (savePSol[i])
2107	oldObj += savePSol[i] * saveCost[i];
2108	if (solution_[i])
2109	newObj += solution_[i] * cost_[i];
2110	bool printIt = false;
2111	if (cost_[i] != saveCost[i])
2112	printIt = true;
2113	if (status_[i] != saveStat[i])
2114	printIt = true;
2115	if (printIt)
2116	printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
2117	i, saveStat[i], saveCost[i], savePSol[i],
2118	status_[i], cost_[i], solution_[i]);
2119	// difference
2120	savePSol[i] = solution_[i] - savePSol[i];
2121	}
2122	printf("exit pivots %d, old obj %g new %g\n",
2123	factorization_->pivots(),
2124	oldObj, newObj);
2125	memset(saveDj, 0, numberRows_ * sizeof(double));
2126	times(1.0, savePSol, saveDj);
2127	double largest = 1.0e-6;
2128	int k = -1;
2129	for (int i = 0; i < numberRows_; i++) {
2130	saveDj[i] -= savePSol[i+numberColumns_];
2131	if (fabs(saveDj[i]) > largest) {
2132	largest = fabs(saveDj[i]);
2133	k = i;
2134	}
2135	}
2136	if (k >= 0)
2137	printf("Not null %d %g\n", k, largest);
2138	}
2139	delete [] savePSol ;
2140	delete [] saveDj ;
2141	delete [] saveCost ;
2142	delete [] saveStat ;
2143	#endif
2144	return returnCode;
2145	}
2146	/* The duals are updated by the given arrays.
2147	Returns number of infeasibilities.
2148	rowArray and columnarray will have flipped
2149	The output vector has movement (row length array) */
2150	int
2151	ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
2152	CoinIndexedVector * columnArray,
2153	CoinIndexedVector * outputArray,
2154	double theta,
2155	double & objectiveChange,
2156	bool fullRecompute)
2157	{
2158
2159	outputArray->clear();
2160
2161
2162	int numberInfeasibilities = 0;
2163	int numberRowInfeasibilities = 0;
2164
2165	// get a tolerance
2166	double tolerance = dualTolerance_;
2167	// we can't really trust infeasibilities if there is dual error
2168	double error = CoinMin(1.0e-2, largestDualError_);
2169	// allow tolerance at least slightly bigger than standard
2170	tolerance = tolerance + error;
2171
2172	double changeObj = 0.0;
2173
2174	// Coding is very similar but we can save a bit by splitting
2175	// Do rows
2176	if (!fullRecompute) {
2177	int i;
2178	double * COIN_RESTRICT reducedCost = djRegion(0);
2179	const double * COIN_RESTRICT lower = lowerRegion(0);
2180	const double * COIN_RESTRICT upper = upperRegion(0);
2181	const double * COIN_RESTRICT cost = costRegion(0);
2182	double * COIN_RESTRICT work;
2183	int number;
2184	int * COIN_RESTRICT which;
2185	const unsigned char * COIN_RESTRICT statusArray = status_ + numberColumns_;
2186	assert(rowArray->packedMode());
2187	work = rowArray->denseVector();
2188	number = rowArray->getNumElements();
2189	which = rowArray->getIndices();
2190	double multiplier[] = { -1.0, 1.0};
2191	for (i = 0; i < number; i++) {
2192	int iSequence = which[i];
2193	double alphaI = work[i];
2194	work[i] = 0.0;
2195	int iStatus = (statusArray[iSequence] & 3) - 1;
2196	if (iStatus) {
2197	double value = reducedCost[iSequence] - theta * alphaI;
2198	reducedCost[iSequence] = value;
2199	double mult = multiplier[iStatus-1];
2200	value *= mult;
2201	if (value < -tolerance) {
2202	// flipping bounds
2203	double movement = mult * (lower[iSequence] - upper[iSequence]);
2204	which[numberInfeasibilities++] = iSequence;
2205	#ifndef NDEBUG
2206	if (fabs(movement) >= 1.0e30)
2207	resetFakeBounds(-1000 - iSequence);
2208	#endif
2209	#ifdef CLP_DEBUG
2210	if ((handler_->logLevel() & 32))
2211	printf("%d %d, new dj %g, alpha %g, movement %g\n",
2212	0, iSequence, value, alphaI, movement);
2213	#endif
2214	changeObj -= movement * cost[iSequence];
2215	outputArray->quickAdd(iSequence, movement);
2216	}
2217	}
2218	}
2219	// Do columns
2220	reducedCost = djRegion(1);
2221	lower = lowerRegion(1);
2222	upper = upperRegion(1);
2223	cost = costRegion(1);
2224	// set number of infeasibilities in row array
2225	numberRowInfeasibilities = numberInfeasibilities;
2226	rowArray->setNumElements(numberInfeasibilities);
2227	numberInfeasibilities = 0;
2228	work = columnArray->denseVector();
2229	number = columnArray->getNumElements();
2230	which = columnArray->getIndices();
2231	if ((moreSpecialOptions_ & 8) != 0) {
2232	const unsigned char * COIN_RESTRICT statusArray = status_;
2233	for (i = 0; i < number; i++) {
2234	int iSequence = which[i];
2235	double alphaI = work[i];
2236	work[i] = 0.0;
2237
2238	int iStatus = (statusArray[iSequence] & 3) - 1;
2239	if (iStatus) {
2240	double value = reducedCost[iSequence] - theta * alphaI;
2241	reducedCost[iSequence] = value;
2242	double mult = multiplier[iStatus-1];
2243	value *= mult;
2244	if (value < -tolerance) {
2245	// flipping bounds
2246	double movement = mult * (upper[iSequence] - lower[iSequence]);
2247	which[numberInfeasibilities++] = iSequence;
2248	#ifndef NDEBUG
2249	if (fabs(movement) >= 1.0e30)
2250	resetFakeBounds(-1000 - iSequence);
2251	#endif
2252	#ifdef CLP_DEBUG
2253	if ((handler_->logLevel() & 32))
2254	printf("%d %d, new dj %g, alpha %g, movement %g\n",
2255	1, iSequence, value, alphaI, movement);
2256	#endif
2257	changeObj += movement * cost[iSequence];
2258	matrix_->add(this, outputArray, iSequence, movement);
2259	}
2260	}
2261	}
2262	} else {
2263	for (i = 0; i < number; i++) {
2264	int iSequence = which[i];
2265	double alphaI = work[i];
2266	work[i] = 0.0;
2267
2268	Status status = getStatus(iSequence);
2269	if (status == atLowerBound) {
2270	double value = reducedCost[iSequence] - theta * alphaI;
2271	reducedCost[iSequence] = value;
2272	double movement = 0.0;
2273
2274	if (value < -tolerance) {
2275	// to upper bound
2276	which[numberInfeasibilities++] = iSequence;
2277	movement = upper[iSequence] - lower[iSequence];
2278	#ifndef NDEBUG
2279	if (fabs(movement) >= 1.0e30)
2280	resetFakeBounds(-1000 - iSequence);
2281	#endif
2282	#ifdef CLP_DEBUG
2283	if ((handler_->logLevel() & 32))
2284	printf("%d %d, new dj %g, alpha %g, movement %g\n",
2285	1, iSequence, value, alphaI, movement);
2286	#endif
2287	changeObj += movement * cost[iSequence];
2288	matrix_->add(this, outputArray, iSequence, movement);
2289	}
2290	} else if (status == atUpperBound) {
2291	double value = reducedCost[iSequence] - theta * alphaI;
2292	reducedCost[iSequence] = value;
2293	double movement = 0.0;
2294
2295	if (value > tolerance) {
2296	// to lower bound (if swap)
2297	which[numberInfeasibilities++] = iSequence;
2298	movement = lower[iSequence] - upper[iSequence];
2299	#ifndef NDEBUG
2300	if (fabs(movement) >= 1.0e30)
2301	resetFakeBounds(-1000 - iSequence);
2302	#endif
2303	#ifdef CLP_DEBUG
2304	if ((handler_->logLevel() & 32))
2305	printf("%d %d, new dj %g, alpha %g, movement %g\n",
2306	1, iSequence, value, alphaI, movement);
2307	#endif
2308	changeObj += movement * cost[iSequence];
2309	matrix_->add(this, outputArray, iSequence, movement);
2310	}
2311	} else if (status == isFree) {
2312	double value = reducedCost[iSequence] - theta * alphaI;
2313	reducedCost[iSequence] = value;
2314	}
2315	}
2316	}
2317	} else {
2318	double * COIN_RESTRICT solution = solutionRegion(0);
2319	double * COIN_RESTRICT reducedCost = djRegion(0);
2320	const double * COIN_RESTRICT lower = lowerRegion(0);
2321	const double * COIN_RESTRICT upper = upperRegion(0);
2322	const double * COIN_RESTRICT cost = costRegion(0);
2323	int * COIN_RESTRICT which;
2324	which = rowArray->getIndices();
2325	int iSequence;
2326	for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2327	double value = reducedCost[iSequence];
2328
2329	Status status = getStatus(iSequence + numberColumns_);
2330	// more likely to be at upper bound ?
2331	if (status == atUpperBound) {
2332	double movement = 0.0;
2333	//#define NO_SWAP7
2334	if (value > tolerance) {
2335	// to lower bound (if swap)
2336	// put back alpha
2337	which[numberInfeasibilities++] = iSequence;
2338	movement = lower[iSequence] - upper[iSequence];
2339	changeObj += movement * cost[iSequence];
2340	outputArray->quickAdd(iSequence, -movement);
2341	#ifndef NDEBUG
2342	if (fabs(movement) >= 1.0e30)
2343	resetFakeBounds(-1000 - iSequence);
2344	#endif
2345	#ifndef NO_SWAP7
2346	} else if (value > -tolerance) {
2347	// at correct bound but may swap
2348	FakeBound bound = getFakeBound(iSequence + numberColumns_);
2349	if (bound == ClpSimplexDual::upperFake) {
2350	movement = lower[iSequence] - upper[iSequence];
2351	#ifndef NDEBUG
2352	if (fabs(movement) >= 1.0e30)
2353	resetFakeBounds(-1000 - iSequence);
2354	#endif
2355	setStatus(iSequence + numberColumns_, atLowerBound);
2356	solution[iSequence] = lower[iSequence];
2357	changeObj += movement * cost[iSequence];
2358	//numberFake_--;
2359	//setFakeBound(iSequence+numberColumns_,noFake);
2360	}
2361	#endif
2362	}
2363	} else if (status == atLowerBound) {
2364	double movement = 0.0;
2365
2366	if (value < -tolerance) {
2367	// to upper bound
2368	// put back alpha
2369	which[numberInfeasibilities++] = iSequence;
2370	movement = upper[iSequence] - lower[iSequence];
2371	#ifndef NDEBUG
2372	if (fabs(movement) >= 1.0e30)
2373	resetFakeBounds(-1000 - iSequence);
2374	#endif
2375	changeObj += movement * cost[iSequence];
2376	outputArray->quickAdd(iSequence, -movement);
2377	#ifndef NO_SWAP7
2378	} else if (value < tolerance) {
2379	// at correct bound but may swap
2380	FakeBound bound = getFakeBound(iSequence + numberColumns_);
2381	if (bound == ClpSimplexDual::lowerFake) {
2382	movement = upper[iSequence] - lower[iSequence];
2383	#ifndef NDEBUG
2384	if (fabs(movement) >= 1.0e30)
2385	resetFakeBounds(-1000 - iSequence);
2386	#endif
2387	setStatus(iSequence + numberColumns_, atUpperBound);
2388	solution[iSequence] = upper[iSequence];
2389	changeObj += movement * cost[iSequence];
2390	//numberFake_--;
2391	//setFakeBound(iSequence+numberColumns_,noFake);
2392	}
2393	#endif
2394	}
2395	}
2396	}
2397	// Do columns
2398	solution = solutionRegion(1);
2399	reducedCost = djRegion(1);
2400	lower = lowerRegion(1);
2401	upper = upperRegion(1);
2402	cost = costRegion(1);
2403	// set number of infeasibilities in row array
2404	numberRowInfeasibilities = numberInfeasibilities;
2405	rowArray->setNumElements(numberInfeasibilities);
2406	numberInfeasibilities = 0;
2407	which = columnArray->getIndices();
2408	for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2409	double value = reducedCost[iSequence];
2410
2411	Status status = getStatus(iSequence);
2412	if (status == atLowerBound) {
2413	double movement = 0.0;
2414
2415	if (value < -tolerance) {
2416	// to upper bound
2417	// put back alpha
2418	which[numberInfeasibilities++] = iSequence;
2419	movement = upper[iSequence] - lower[iSequence];
2420	#ifndef NDEBUG
2421	if (fabs(movement) >= 1.0e30)
2422	resetFakeBounds(-1000 - iSequence);
2423	#endif
2424	changeObj += movement * cost[iSequence];
2425	matrix_->add(this, outputArray, iSequence, movement);
2426	#ifndef NO_SWAP7
2427	} else if (value < tolerance) {
2428	// at correct bound but may swap
2429	FakeBound bound = getFakeBound(iSequence);
2430	if (bound == ClpSimplexDual::lowerFake) {
2431	movement = upper[iSequence] - lower[iSequence];
2432	#ifndef NDEBUG
2433	if (fabs(movement) >= 1.0e30)
2434	resetFakeBounds(-1000 - iSequence);
2435	#endif
2436	setStatus(iSequence, atUpperBound);
2437	solution[iSequence] = upper[iSequence];
2438	changeObj += movement * cost[iSequence];
2439	//numberFake_--;
2440	//setFakeBound(iSequence,noFake);
2441	}
2442	#endif
2443	}
2444	} else if (status == atUpperBound) {
2445	double movement = 0.0;
2446
2447	if (value > tolerance) {
2448	// to lower bound (if swap)
2449	// put back alpha
2450	which[numberInfeasibilities++] = iSequence;
2451	movement = lower[iSequence] - upper[iSequence];
2452	#ifndef NDEBUG
2453	if (fabs(movement) >= 1.0e30)
2454	resetFakeBounds(-1000 - iSequence);
2455	#endif
2456	changeObj += movement * cost[iSequence];
2457	matrix_->add(this, outputArray, iSequence, movement);
2458	#ifndef NO_SWAP7
2459	} else if (value > -tolerance) {
2460	// at correct bound but may swap
2461	FakeBound bound = getFakeBound(iSequence);
2462	if (bound == ClpSimplexDual::upperFake) {
2463	movement = lower[iSequence] - upper[iSequence];
2464	#ifndef NDEBUG
2465	if (fabs(movement) >= 1.0e30)
2466	resetFakeBounds(-1000 - iSequence);
2467	#endif
2468	setStatus(iSequence, atLowerBound);
2469	solution[iSequence] = lower[iSequence];
2470	changeObj += movement * cost[iSequence];
2471	//numberFake_--;
2472	//setFakeBound(iSequence,noFake);
2473	}
2474	#endif
2475	}
2476	}
2477	}
2478	}
2479
2480	#ifdef CLP_DEBUG
2481	if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0)
2482	printf("%d fake after full update\n", numberFake_);
2483	#endif
2484	// set number of infeasibilities
2485	columnArray->setNumElements(numberInfeasibilities);
2486	numberInfeasibilities += numberRowInfeasibilities;
2487	if (fullRecompute) {
2488	// do actual flips
2489	flipBounds(rowArray, columnArray);
2490	}
2491	objectiveChange += changeObj;
2492	return numberInfeasibilities;
2493	}
2494	void
2495	ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
2496	CoinIndexedVector * columnArray,
2497	double theta)
2498	{
2499
2500	// use a tighter tolerance except for all being okay
2501	double tolerance = dualTolerance_;
2502
2503	// Coding is very similar but we can save a bit by splitting
2504	// Do rows
2505	{
2506	int i;
2507	double * reducedCost = djRegion(0);
2508	double * work;
2509	int number;
2510	int * which;
2511	work = rowArray->denseVector();
2512	number = rowArray->getNumElements();
2513	which = rowArray->getIndices();
2514	for (i = 0; i < number; i++) {
2515	int iSequence = which[i];
2516	double alphaI = work[i];
2517	double value = reducedCost[iSequence] - theta * alphaI;
2518	work[i] = 0.0;
2519	reducedCost[iSequence] = value;
2520
2521	Status status = getStatus(iSequence + numberColumns_);
2522	// more likely to be at upper bound ?
2523	if (status == atUpperBound) {
2524
2525	if (value > tolerance)
2526	reducedCost[iSequence] = 0.0;
2527	} else if (status == atLowerBound) {
2528
2529	if (value < -tolerance) {
2530	reducedCost[iSequence] = 0.0;
2531	}
2532	}
2533	}
2534	}
2535	rowArray->setNumElements(0);
2536
2537	// Do columns
2538	{
2539	int i;
2540	double * reducedCost = djRegion(1);
2541	double * work;
2542	int number;
2543	int * which;
2544	work = columnArray->denseVector();
2545	number = columnArray->getNumElements();
2546	which = columnArray->getIndices();
2547
2548	for (i = 0; i < number; i++) {
2549	int iSequence = which[i];
2550	double alphaI = work[i];
2551	double value = reducedCost[iSequence] - theta * alphaI;
2552	work[i] = 0.0;
2553	reducedCost[iSequence] = value;
2554
2555	Status status = getStatus(iSequence);
2556	if (status == atLowerBound) {
2557	if (value < -tolerance)
2558	reducedCost[iSequence] = 0.0;
2559	} else if (status == atUpperBound) {
2560	if (value > tolerance)
2561	reducedCost[iSequence] = 0.0;
2562	}
2563	}
2564	}
2565	columnArray->setNumElements(0);
2566	}
2567	/*
2568	Chooses dual pivot row
2569	Would be faster with separate region to scan
2570	and will have this (with square of infeasibility) when steepest
2571	For easy problems we can just choose one of the first rows we look at
2572	*/
2573	void
2574	ClpSimplexDual::dualRow(int alreadyChosen)
2575	{
2576	// get pivot row using whichever method it is
2577	int chosenRow = -1;
2578	#ifdef FORCE_FOLLOW
2579	bool forceThis = false;
2580	if (!fpFollow && strlen(forceFile)) {
2581	fpFollow = fopen(forceFile, "r");
2582	assert (fpFollow);
2583	}
2584	if (fpFollow) {
2585	if (numberIterations_ <= force_iteration) {
2586	// read to next Clp0102
2587	char temp[300];
2588	while (fgets(temp, 250, fpFollow)) {
2589	if (strncmp(temp, "Clp0102", 7))
2590	continue;
2591	char cin, cout;
2592	sscanf(temp + 9, "%d%f%s%c%c%d%s%*c%c%d",
2593	&force_iteration, &cin, &force_in, &cout, &force_out);
2594	if (cin == 'R')
2595	force_in += numberColumns_;
2596	if (cout == 'R')
2597	force_out += numberColumns_;
2598	forceThis = true;
2599	assert (numberIterations_ == force_iteration - 1);
2600	printf("Iteration %d will force %d out and %d in\n",
2601	force_iteration, force_out, force_in);
2602	alreadyChosen = force_out;
2603	break;
2604	}
2605	} else {
2606	// use old
2607	forceThis = true;
2608	}
2609	if (!forceThis) {
2610	fclose(fpFollow);
2611	fpFollow = NULL;
2612	forceFile = "";
2613	}
2614	}
2615	#endif
2616	double freeAlpha = 0.0;
2617	if (alreadyChosen < 0) {
2618	// first see if any free variables and put them in basis
2619	int nextFree = nextSuperBasic();
2620	//nextFree=-1; //off
2621	if (nextFree >= 0) {
2622	// unpack vector and find a good pivot
2623	unpack(rowArray_[1], nextFree);
2624	factorization_->updateColumn(rowArray_[2], rowArray_[1]);
2625
2626	double * work = rowArray_[1]->denseVector();
2627	int number = rowArray_[1]->getNumElements();
2628	int * which = rowArray_[1]->getIndices();
2629	double bestFeasibleAlpha = 0.0;
2630	int bestFeasibleRow = -1;
2631	double bestInfeasibleAlpha = 0.0;
2632	int bestInfeasibleRow = -1;
2633	int i;
2634
2635	for (i = 0; i < number; i++) {
2636	int iRow = which[i];
2637	double alpha = fabs(work[iRow]);
2638	if (alpha > 1.0e-3) {
2639	int iSequence = pivotVariable_[iRow];
2640	double value = solution_[iSequence];
2641	double lower = lower_[iSequence];
2642	double upper = upper_[iSequence];
2643	double infeasibility = 0.0;
2644	if (value > upper)
2645	infeasibility = value - upper;
2646	else if (value < lower)
2647	infeasibility = lower - value;
2648	if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
2649	if (!flagged(iSequence)) {
2650	bestInfeasibleAlpha = infeasibility * alpha;
2651	bestInfeasibleRow = iRow;
2652	}
2653	}
2654	if (alpha > bestFeasibleAlpha && (lower > -1.0e20 \|\| upper < 1.0e20)) {
2655	bestFeasibleAlpha = alpha;
2656	bestFeasibleRow = iRow;
2657	}
2658	}
2659	}
2660	if (bestInfeasibleRow >= 0)
2661	chosenRow = bestInfeasibleRow;
2662	else if (bestFeasibleAlpha > 1.0e-2)
2663	chosenRow = bestFeasibleRow;
2664	if (chosenRow >= 0) {
2665	pivotRow_ = chosenRow;
2666	freeAlpha = work[chosenRow];
2667	}
2668	rowArray_[1]->clear();
2669	}
2670	} else {
2671	// in values pass
2672	chosenRow = alreadyChosen;
2673	#ifdef FORCE_FOLLOW
2674	if(forceThis) {
2675	alreadyChosen = -1;
2676	chosenRow = -1;
2677	for (int i = 0; i < numberRows_; i++) {
2678	if (pivotVariable_[i] == force_out) {
2679	chosenRow = i;
2680	break;
2681	}
2682	}
2683	assert (chosenRow >= 0);
2684	}
2685	#endif
2686	pivotRow_ = chosenRow;
2687	}
2688	if (chosenRow < 0)
2689	pivotRow_ = dualRowPivot_->pivotRow();
2690
2691	if (pivotRow_ >= 0) {
2692	sequenceOut_ = pivotVariable_[pivotRow_];
2693	valueOut_ = solution_[sequenceOut_];
2694	lowerOut_ = lower_[sequenceOut_];
2695	upperOut_ = upper_[sequenceOut_];
2696	if (alreadyChosen < 0) {
2697	// if we have problems we could try other way and hope we get a
2698	// zero pivot?
2699	if (valueOut_ > upperOut_) {
2700	directionOut_ = -1;
2701	dualOut_ = valueOut_ - upperOut_;
2702	} else if (valueOut_ < lowerOut_) {
2703	directionOut_ = 1;
2704	dualOut_ = lowerOut_ - valueOut_;
2705	} else {
2706	#if 1
2707	// odd (could be free) - it's feasible - go to nearest
2708	if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
2709	directionOut_ = 1;
2710	dualOut_ = lowerOut_ - valueOut_;
2711	} else {
2712	directionOut_ = -1;
2713	dualOut_ = valueOut_ - upperOut_;
2714	}
2715	#else
2716	// odd (could be free) - it's feasible - improve obj
2717	printf("direction from alpha of %g is %d\n",
2718	freeAlpha, freeAlpha > 0.0 ? 1 : -1);
2719	if (valueOut_ - lowerOut_ > 1.0e20)
2720	freeAlpha = 1.0;
2721	else if(upperOut_ - valueOut_ > 1.0e20)
2722	freeAlpha = -1.0;
2723	//if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2724	if (freeAlpha < 0.0) {
2725	directionOut_ = 1;
2726	dualOut_ = lowerOut_ - valueOut_;
2727	} else {
2728	directionOut_ = -1;
2729	dualOut_ = valueOut_ - upperOut_;
2730	}
2731	printf("direction taken %d - bounds %g %g %g\n",
2732	directionOut_, lowerOut_, valueOut_, upperOut_);
2733	#endif
2734	}
2735	#ifdef CLP_DEBUG
2736	assert(dualOut_ >= 0.0);
2737	#endif
2738	} else {
2739	// in values pass so just use sign of dj
2740	// We don't want to go through any barriers so set dualOut low
2741	// free variables will never be here
2742	dualOut_ = 1.0e-6;
2743	if (dj_[sequenceOut_] > 0.0) {
2744	// this will give a -1 in pivot row (as slacks are -1.0)
2745	directionOut_ = 1;
2746	} else {
2747	directionOut_ = -1;
2748	}
2749	}
2750	}
2751	return ;
2752	}
2753	// Checks if any fake bounds active - if so returns number and modifies
2754	// dualBound_ and everything.
2755	// Free variables will be left as free
2756	// Returns number of bounds changed if >=0
2757	// Returns -1 if not initialize and no effect
2758	// Fills in changeVector which can be used to see if unbounded
2759	// and cost of change vector
2760	int
2761	ClpSimplexDual::changeBounds(int initialize,
2762	CoinIndexedVector * outputArray,
2763	double & changeCost)
2764	{
2765	numberFake_ = 0;
2766	if (!initialize) {
2767	int numberInfeasibilities;
2768	double newBound;
2769	newBound = 5.0 * dualBound_;
2770	numberInfeasibilities = 0;
2771	changeCost = 0.0;
2772	// put back original bounds and then check
2773	createRim1(false);
2774	int iSequence;
2775	// bounds will get bigger - just look at ones at bounds
2776	for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2777	double lowerValue = lower_[iSequence];
2778	double upperValue = upper_[iSequence];
2779	double value = solution_[iSequence];
2780	setFakeBound(iSequence, ClpSimplexDual::noFake);
2781	switch(getStatus(iSequence)) {
2782
2783	case basic:
2784	case ClpSimplex::isFixed:
2785	break;
2786	case isFree:
2787	case superBasic:
2788	break;
2789	case atUpperBound:
2790	if (fabs(value - upperValue) > primalTolerance_)
2791	numberInfeasibilities++;
2792	break;
2793	case atLowerBound:
2794	if (fabs(value - lowerValue) > primalTolerance_)
2795	numberInfeasibilities++;
2796	break;
2797	}
2798	}
2799	// If dual infeasible then carry on
2800	if (numberInfeasibilities) {
2801	handler_->message(CLP_DUAL_CHECKB, messages_)
2802	<< newBound
2803	<< CoinMessageEol;
2804	int iSequence;
2805	for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2806	double lowerValue = lower_[iSequence];
2807	double upperValue = upper_[iSequence];
2808	double newLowerValue;
2809	double newUpperValue;
2810	Status status = getStatus(iSequence);
2811	if (status == atUpperBound \|\|
2812	status == atLowerBound) {
2813	double value = solution_[iSequence];
2814	if (value - lowerValue <= upperValue - value) {
2815	newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound);
2816	newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
2817	} else {
2818	newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound);
2819	newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
2820	}
2821	lower_[iSequence] = newLowerValue;
2822	upper_[iSequence] = newUpperValue;
2823	if (newLowerValue > lowerValue) {
2824	if (newUpperValue < upperValue) {
2825	setFakeBound(iSequence, ClpSimplexDual::bothFake);
2826	#ifdef CLP_INVESTIGATE
2827	abort(); // No idea what should happen here - I have never got here
2828	#endif
2829	numberFake_++;
2830	} else {
2831	setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2832	numberFake_++;
2833	}
2834	} else {
2835	if (newUpperValue < upperValue) {
2836	setFakeBound(iSequence, ClpSimplexDual::upperFake);
2837	numberFake_++;
2838	}
2839	}
2840	if (status == atUpperBound)
2841	solution_[iSequence] = newUpperValue;
2842	else
2843	solution_[iSequence] = newLowerValue;
2844	double movement = solution_[iSequence] - value;
2845	if (movement && outputArray) {
2846	if (iSequence >= numberColumns_) {
2847	outputArray->quickAdd(iSequence, -movement);
2848	changeCost += movement * cost_[iSequence];
2849	} else {
2850	matrix_->add(this, outputArray, iSequence, movement);
2851	changeCost += movement * cost_[iSequence];
2852	}
2853	}
2854	}
2855	}
2856	dualBound_ = newBound;
2857	} else {
2858	numberInfeasibilities = -1;
2859	}
2860	return numberInfeasibilities;
2861	} else if (initialize == 1 \|\| initialize == 3) {
2862	int iSequence;
2863	if (initialize == 3) {
2864	for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2865	setFakeBound(iSequence, ClpSimplexDual::noFake);
2866	}
2867	}
2868	double testBound = 0.999999 * dualBound_;
2869	for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2870	Status status = getStatus(iSequence);
2871	if (status == atUpperBound \|\|
2872	status == atLowerBound) {
2873	double lowerValue = lower_[iSequence];
2874	double upperValue = upper_[iSequence];
2875	double value = solution_[iSequence];
2876	if (lowerValue > -largeValue_ \|\| upperValue < largeValue_) {
2877	if (true \|\| lowerValue - value > -0.5 * dualBound_ \|\|
2878	upperValue - value < 0.5 * dualBound_) {
2879	if (fabs(lowerValue - value) <= fabs(upperValue - value)) {
2880	if (upperValue > lowerValue + testBound) {
2881	if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2882	numberFake_++;
2883	upper_[iSequence] = lowerValue + dualBound_;
2884	setFakeBound(iSequence, ClpSimplexDual::upperFake);
2885	}
2886	} else {
2887	if (lowerValue < upperValue - testBound) {
2888	if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2889	numberFake_++;
2890	lower_[iSequence] = upperValue - dualBound_;
2891	setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2892	}
2893	}
2894	} else {
2895	if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2896	numberFake_++;
2897	lower_[iSequence] = -0.5 * dualBound_;
2898	upper_[iSequence] = 0.5 * dualBound_;
2899	setFakeBound(iSequence, ClpSimplexDual::bothFake);
2900	abort();
2901	}
2902	if (status == atUpperBound)
2903	solution_[iSequence] = upper_[iSequence];
2904	else
2905	solution_[iSequence] = lower_[iSequence];
2906	} else {
2907	// set non basic free variables to fake bounds
2908	// I don't think we should ever get here
2909	CoinAssert(!("should not be here"));
2910	lower_[iSequence] = -0.5 * dualBound_;
2911	upper_[iSequence] = 0.5 * dualBound_;
2912	setFakeBound(iSequence, ClpSimplexDual::bothFake);
2913	numberFake_++;
2914	setStatus(iSequence, atUpperBound);
2915	solution_[iSequence] = 0.5 * dualBound_;
2916	}
2917	} else if (status == basic) {
2918	// make sure not at fake bound and bounds correct
2919	setFakeBound(iSequence, ClpSimplexDual::noFake);
2920	double gap = upper_[iSequence] - lower_[iSequence];
2921	if (gap > 0.5 * dualBound_ && gap < 2.0 * dualBound_) {
2922	if (iSequence < numberColumns_) {
2923	if (columnScale_) {
2924	double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2925	// lower
2926	double value = columnLower_[iSequence];
2927	if (value > -1.0e30) {
2928	value *= multiplier;
2929	}
2930	lower_[iSequence] = value;
2931	// upper
2932	value = columnUpper_[iSequence];
2933	if (value < 1.0e30) {
2934	value *= multiplier;
2935	}
2936	upper_[iSequence] = value;
2937	} else {
2938	lower_[iSequence] = columnLower_[iSequence];;
2939	upper_[iSequence] = columnUpper_[iSequence];;
2940	}
2941	} else {
2942	int iRow = iSequence - numberColumns_;
2943	if (rowScale_) {
2944	// lower
2945	double multiplier = rhsScale_ * rowScale_[iRow];
2946	double value = rowLower_[iRow];
2947	if (value > -1.0e30) {
2948	value *= multiplier;
2949	}
2950	lower_[iSequence] = value;
2951	// upper
2952	value = rowUpper_[iRow];
2953	if (value < 1.0e30) {
2954	value *= multiplier;
2955	}
2956	upper_[iSequence] = value;
2957	} else {
2958	lower_[iSequence] = rowLower_[iRow];;
2959	upper_[iSequence] = rowUpper_[iRow];;
2960	}
2961	}
2962	}
2963	}
2964	}
2965
2966	return 1;
2967	} else {
2968	// just reset changed ones
2969	if (columnScale_) {
2970	int iSequence;
2971	for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2972	FakeBound fakeStatus = getFakeBound(iSequence);
2973	if (fakeStatus != noFake) {
2974	if ((static_cast<int> (fakeStatus) & 1) != 0) {
2975	// lower
2976	double value = columnLower_[iSequence];
2977	if (value > -1.0e30) {
2978	double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2979	value *= multiplier;
2980	}
2981	columnLowerWork_[iSequence] = value;
2982	}
2983	if ((static_cast<int> (fakeStatus) & 2) != 0) {
2984	// upper
2985	double value = columnUpper_[iSequence];
2986	if (value < 1.0e30) {
2987	double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2988	value *= multiplier;
2989	}
2990	columnUpperWork_[iSequence] = value;
2991	}
2992	}
2993	}
2994	for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2995	FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
2996	if (fakeStatus != noFake) {
2997	if ((static_cast<int> (fakeStatus) & 1) != 0) {
2998	// lower
2999	double value = rowLower_[iSequence];
3000	if (value > -1.0e30) {
3001	double multiplier = rhsScale_ * rowScale_[iSequence];
3002	value *= multiplier;
3003	}
3004	rowLowerWork_[iSequence] = value;
3005	}
3006	if ((static_cast<int> (fakeStatus) & 2) != 0) {
3007	// upper
3008	double value = rowUpper_[iSequence];
3009	if (value < 1.0e30) {
3010	double multiplier = rhsScale_ * rowScale_[iSequence];
3011	value *= multiplier;
3012	}
3013	rowUpperWork_[iSequence] = value;
3014	}
3015	}
3016	}
3017	} else {
3018	int iSequence;
3019	for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3020	FakeBound fakeStatus = getFakeBound(iSequence);
3021	if ((static_cast<int> (fakeStatus) & 1) != 0) {
3022	// lower
3023	columnLowerWork_[iSequence] = columnLower_[iSequence];
3024	}
3025	if ((static_cast<int> (fakeStatus) & 2) != 0) {
3026	// upper
3027	columnUpperWork_[iSequence] = columnUpper_[iSequence];
3028	}
3029	}
3030	for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3031	FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3032	if ((static_cast<int> (fakeStatus) & 1) != 0) {
3033	// lower
3034	rowLowerWork_[iSequence] = rowLower_[iSequence];
3035	}
3036	if ((static_cast<int> (fakeStatus) & 2) != 0) {
3037	// upper
3038	rowUpperWork_[iSequence] = rowUpper_[iSequence];
3039	}
3040	}
3041	}
3042	return 0;
3043	}
3044	}
3045	int
3046	ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
3047	const CoinIndexedVector * columnArray,
3048	CoinIndexedVector * spareArray,
3049	double acceptablePivot,
3050	double & upperReturn, double &bestReturn, double & badFree)
3051	{
3052	// do first pass to get possibles
3053	double * spare = spareArray->denseVector();
3054	int * index = spareArray->getIndices();
3055	const double * work;
3056	int number;
3057	const int * which;
3058	const double * reducedCost;
3059	// We can also see if infeasible or pivoting on free
3060	double tentativeTheta = 1.0e15;
3061	double upperTheta = 1.0e31;
3062	double freePivot = acceptablePivot;
3063	double bestPossible = 0.0;
3064	int numberRemaining = 0;
3065	int i;
3066	badFree = 0.0;
3067	if ((moreSpecialOptions_ & 8) != 0) {
3068	// No free or super basic
3069	double multiplier[] = { -1.0, 1.0};
3070	double dualT = - dualTolerance_;
3071	for (int iSection = 0; iSection < 2; iSection++) {
3072
3073	int addSequence;
3074	unsigned char * statusArray;
3075	if (!iSection) {
3076	work = rowArray->denseVector();
3077	number = rowArray->getNumElements();
3078	which = rowArray->getIndices();
3079	reducedCost = rowReducedCost_;
3080	addSequence = numberColumns_;
3081	statusArray = status_ + numberColumns_;
3082	} else {
3083	work = columnArray->denseVector();
3084	number = columnArray->getNumElements();
3085	which = columnArray->getIndices();
3086	reducedCost = reducedCostWork_;
3087	addSequence = 0;
3088	statusArray = status_;
3089	}
3090
3091	for (i = 0; i < number; i++) {
3092	int iSequence = which[i];
3093	double alpha;
3094	double oldValue;
3095	double value;
3096
3097	assert (getStatus(iSequence + addSequence) != isFree
3098	&& getStatus(iSequence + addSequence) != superBasic);
3099	int iStatus = (statusArray[iSequence] & 3) - 1;
3100	if (iStatus) {
3101	double mult = multiplier[iStatus-1];
3102	alpha = work[i] * mult;
3103	if (alpha > 0.0) {
3104	oldValue = reducedCost[iSequence] * mult;
3105	value = oldValue - tentativeTheta * alpha;
3106	if (value < dualT) {
3107	bestPossible = CoinMax(bestPossible, alpha);
3108	value = oldValue - upperTheta * alpha;
3109	if (value < dualT && alpha >= acceptablePivot) {
3110	upperTheta = (oldValue - dualT) / alpha;
3111	//tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3112	}
3113	// add to list
3114	spare[numberRemaining] = alpha * mult;
3115	index[numberRemaining++] = iSequence + addSequence;
3116	}
3117	}
3118	}
3119	}
3120	}
3121	} else {
3122	// some free or super basic
3123	for (int iSection = 0; iSection < 2; iSection++) {
3124
3125	int addSequence;
3126
3127	if (!iSection) {
3128	work = rowArray->denseVector();
3129	number = rowArray->getNumElements();
3130	which = rowArray->getIndices();
3131	reducedCost = rowReducedCost_;
3132	addSequence = numberColumns_;
3133	} else {
3134	work = columnArray->denseVector();
3135	number = columnArray->getNumElements();
3136	which = columnArray->getIndices();
3137	reducedCost = reducedCostWork_;
3138	addSequence = 0;
3139	}
3140
3141	for (i = 0; i < number; i++) {
3142	int iSequence = which[i];
3143	double alpha;
3144	double oldValue;
3145	double value;
3146	bool keep;
3147
3148	switch(getStatus(iSequence + addSequence)) {
3149
3150	case basic:
3151	case ClpSimplex::isFixed:
3152	break;
3153	case isFree:
3154	case superBasic:
3155	alpha = work[i];
3156	bestPossible = CoinMax(bestPossible, fabs(alpha));
3157	oldValue = reducedCost[iSequence];
3158	// If free has to be very large - should come in via dualRow
3159	//if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
3160	//break;
3161	if (oldValue > dualTolerance_) {
3162	keep = true;
3163	} else if (oldValue < -dualTolerance_) {
3164	keep = true;
3165	} else {
3166	if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
3167	keep = true;
3168	} else {
3169	keep = false;
3170	badFree = CoinMax(badFree, fabs(alpha));
3171	}
3172	}
3173	if (keep) {
3174	// free - choose largest
3175	if (fabs(alpha) > freePivot) {
3176	freePivot = fabs(alpha);
3177	sequenceIn_ = iSequence + addSequence;
3178	theta_ = oldValue / alpha;
3179	alpha_ = alpha;
3180	}
3181	}
3182	break;
3183	case atUpperBound:
3184	alpha = work[i];
3185	oldValue = reducedCost[iSequence];
3186	value = oldValue - tentativeTheta * alpha;
3187	//assert (oldValue<=dualTolerance_*1.0001);
3188	if (value > dualTolerance_) {
3189	bestPossible = CoinMax(bestPossible, -alpha);
3190	value = oldValue - upperTheta * alpha;
3191	if (value > dualTolerance_ && -alpha >= acceptablePivot) {
3192	upperTheta = (oldValue - dualTolerance_) / alpha;
3193	//tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3194	}
3195	// add to list
3196	spare[numberRemaining] = alpha;
3197	index[numberRemaining++] = iSequence + addSequence;
3198	}
3199	break;
3200	case atLowerBound:
3201	alpha = work[i];
3202	oldValue = reducedCost[iSequence];
3203	value = oldValue - tentativeTheta * alpha;
3204	//assert (oldValue>=-dualTolerance_*1.0001);
3205	if (value < -dualTolerance_) {
3206	bestPossible = CoinMax(bestPossible, alpha);
3207	value = oldValue - upperTheta * alpha;
3208	if (value < -dualTolerance_ && alpha >= acceptablePivot) {
3209	upperTheta = (oldValue + dualTolerance_) / alpha;
3210	//tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3211	}
3212	// add to list
3213	spare[numberRemaining] = alpha;
3214	index[numberRemaining++] = iSequence + addSequence;
3215	}
3216	break;
3217	}
3218	}
3219	}
3220	}
3221	upperReturn = upperTheta;
3222	bestReturn = bestPossible;
3223	return numberRemaining;
3224	}
3225	/*
3226	Row array has row part of pivot row (as duals so sign may be switched)
3227	Column array has column part.
3228	This chooses pivot column.
3229	Spare array will be needed when we start getting clever.
3230	We will check for basic so spare array will never overflow.
3231	If necessary will modify costs
3232	*/
3233	double
3234	ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
3235	CoinIndexedVector * columnArray,
3236	CoinIndexedVector * spareArray,
3237	CoinIndexedVector * spareArray2,
3238	double acceptablePivot,
3239	CoinBigIndex * /dubiousWeights/)
3240	{
3241	int numberPossiblySwapped = 0;
3242	int numberRemaining = 0;
3243
3244	double totalThru = 0.0; // for when variables flip
3245	//double saveAcceptable=acceptablePivot;
3246	//acceptablePivot=1.0e-9;
3247
3248	double bestEverPivot = acceptablePivot;
3249	int lastSequence = -1;
3250	double lastPivot = 0.0;
3251	double upperTheta;
3252	double newTolerance = dualTolerance_;
3253	//newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
3254	// will we need to increase tolerance
3255	bool thisIncrease = false;
3256	// If we think we need to modify costs (not if something from broad sweep)
3257	bool modifyCosts = false;
3258	// Increase in objective due to swapping bounds (may be negative)
3259	double increaseInObjective = 0.0;
3260
3261	// use spareArrays to put ones looked at in
3262	// we are going to flip flop between
3263	int iFlip = 0;
3264	// Possible list of pivots
3265	int interesting[2];
3266	// where possible swapped ones are
3267	int swapped[2];
3268	// for zeroing out arrays after
3269	int marker[2][2];
3270	// pivot elements
3271	double * array[2], * spare, * spare2;
3272	// indices
3273	int * indices[2], * index, * index2;
3274	spareArray2->clear();
3275	array[0] = spareArray->denseVector();
3276	indices[0] = spareArray->getIndices();
3277	spare = array[0];
3278	index = indices[0];
3279	array[1] = spareArray2->denseVector();
3280	indices[1] = spareArray2->getIndices();
3281	int i;
3282
3283	// initialize lists
3284	for (i = 0; i < 2; i++) {
3285	interesting[i] = 0;
3286	swapped[i] = numberColumns_;
3287	marker[i][0] = 0;
3288	marker[i][1] = numberColumns_;
3289	}
3290	/*
3291	First we get a list of possible pivots. We can also see if the
3292	problem looks infeasible or whether we want to pivot in free variable.
3293	This may make objective go backwards but can only happen a finite
3294	number of times and I do want free variables basic.
3295
3296	Then we flip back and forth. At the start of each iteration
3297	interesting[iFlip] should have possible candidates and swapped[iFlip]
3298	will have pivots if we decide to take a previous pivot.
3299	At end of each iteration interesting[1-iFlip] should have
3300	candidates if we go through this theta and swapped[1-iFlip]
3301	pivots if we don't go through.
3302
3303	At first we increase theta and see what happens. We start
3304	theta at a reasonable guess. If in right area then we do bit by bit.
3305
3306	*/
3307
3308	// do first pass to get possibles
3309	upperTheta = 1.0e31;
3310	double bestPossible = 0.0;
3311	double badFree = 0.0;
3312	alpha_ = 0.0;
3313	if (spareIntArray_[0] >= 0) {
3314	numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3315	acceptablePivot, upperTheta, bestPossible, badFree);
3316	} else {
3317	// already done
3318	numberRemaining = spareArray->getNumElements();
3319	spareArray->setNumElements(0);
3320	upperTheta = spareDoubleArray_[0];
3321	bestPossible = spareDoubleArray_[1];
3322	if (spareIntArray_[0] == -1) {
3323	theta_ = spareDoubleArray_[2];
3324	alpha_ = spareDoubleArray_[3];
3325	sequenceIn_ = spareIntArray_[1];
3326	} else {
3327	#if 0
3328	int n = numberRemaining;
3329	double u = upperTheta;
3330	double b = bestPossible;
3331	upperTheta = 1.0e31;
3332	bestPossible = 0.0;
3333	numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3334	acceptablePivot, upperTheta, bestPossible, badFree);
3335	assert (n == numberRemaining);
3336	assert (fabs(b - bestPossible) < 1.0e-7);
3337	assert (fabs(u - upperTheta) < 1.0e-7);
3338	#endif
3339	}
3340	}
3341	// switch off
3342	spareIntArray_[0] = 0;
3343	// We can also see if infeasible or pivoting on free
3344	double tentativeTheta = 1.0e25;
3345	interesting[0] = numberRemaining;
3346	marker[0][0] = numberRemaining;
3347
3348	if (!numberRemaining && sequenceIn_ < 0)
3349	return 0.0; // Looks infeasible
3350
3351	// If sum of bad small pivots too much
3352	#define MORE_CAREFUL
3353	#ifdef MORE_CAREFUL
3354	bool badSumPivots = false;
3355	#endif
3356	if (sequenceIn_ >= 0) {
3357	// free variable - always choose
3358	} else {
3359
3360	theta_ = 1.0e50;
3361	// now flip flop between spare arrays until reasonable theta
3362	tentativeTheta = CoinMax(10.0 * upperTheta, 1.0e-7);
3363
3364	// loops increasing tentative theta until can't go through
3365
3366	while (tentativeTheta < 1.0e22) {
3367	double thruThis = 0.0;
3368
3369	double bestPivot = acceptablePivot;
3370	int bestSequence = -1;
3371
3372	numberPossiblySwapped = numberColumns_;
3373	numberRemaining = 0;
3374
3375	upperTheta = 1.0e50;
3376
3377	spare = array[iFlip];
3378	index = indices[iFlip];
3379	spare2 = array[1-iFlip];
3380	index2 = indices[1-iFlip];
3381
3382	// try 3 different ways
3383	// 1 bias increase by ones with slightly wrong djs
3384	// 2 bias by all
3385	// 3 bias by all - tolerance
3386	#define TRYBIAS 3
3387
3388
3389	double increaseInThis = 0.0; //objective increase in this loop
3390
3391	for (i = 0; i < interesting[iFlip]; i++) {
3392	int iSequence = index[i];
3393	double alpha = spare[i];
3394	double oldValue = dj_[iSequence];
3395	double value = oldValue - tentativeTheta * alpha;
3396
3397	if (alpha < 0.0) {
3398	//at upper bound
3399	if (value > newTolerance) {
3400	double range = upper_[iSequence] - lower_[iSequence];
3401	thruThis -= range * alpha;
3402	#if TRYBIAS==1
3403	if (oldValue > 0.0)
3404	increaseInThis -= oldValue * range;
3405	#elif TRYBIAS==2
3406	increaseInThis -= oldValue * range;
3407	#else
3408	increaseInThis -= (oldValue + dualTolerance_) * range;
3409	#endif
3410	// goes on swapped list (also means candidates if too many)
3411	spare2[--numberPossiblySwapped] = alpha;
3412	index2[numberPossiblySwapped] = iSequence;
3413	if (fabs(alpha) > bestPivot) {
3414	bestPivot = fabs(alpha);
3415	bestSequence = numberPossiblySwapped;
3416	}
3417	} else {
3418	value = oldValue - upperTheta * alpha;
3419	if (value > newTolerance && -alpha >= acceptablePivot)
3420	upperTheta = (oldValue - newTolerance) / alpha;
3421	spare2[numberRemaining] = alpha;
3422	index2[numberRemaining++] = iSequence;
3423	}
3424	} else {
3425	// at lower bound
3426	if (value < -newTolerance) {
3427	double range = upper_[iSequence] - lower_[iSequence];
3428	thruThis += range * alpha;
3429	//?? is this correct - and should we look at good ones
3430	#if TRYBIAS==1
3431	if (oldValue < 0.0)
3432	increaseInThis += oldValue * range;
3433	#elif TRYBIAS==2
3434	increaseInThis += oldValue * range;
3435	#else
3436	increaseInThis += (oldValue - dualTolerance_) * range;
3437	#endif
3438	// goes on swapped list (also means candidates if too many)
3439	spare2[--numberPossiblySwapped] = alpha;
3440	index2[numberPossiblySwapped] = iSequence;
3441	if (fabs(alpha) > bestPivot) {
3442	bestPivot = fabs(alpha);
3443	bestSequence = numberPossiblySwapped;
3444	}
3445	} else {
3446	value = oldValue - upperTheta * alpha;
3447	if (value < -newTolerance && alpha >= acceptablePivot)
3448	upperTheta = (oldValue + newTolerance) / alpha;
3449	spare2[numberRemaining] = alpha;
3450	index2[numberRemaining++] = iSequence;
3451	}
3452	}
3453	}
3454	swapped[1-iFlip] = numberPossiblySwapped;
3455	interesting[1-iFlip] = numberRemaining;
3456	marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3457	marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3458
3459	if (totalThru + thruThis >= fabs(dualOut_) \|\|
3460	increaseInObjective + increaseInThis < 0.0) {
3461	// We should be pivoting in this batch
3462	// so compress down to this lot
3463	numberRemaining = 0;
3464	for (i = numberColumns_ - 1; i >= swapped[1-iFlip]; i--) {
3465	spare[numberRemaining] = spare2[i];
3466	index[numberRemaining++] = index2[i];
3467	}
3468	interesting[iFlip] = numberRemaining;
3469	int iTry;
3470	#define MAXTRY 100
3471	// first get ratio with tolerance
3472	for (iTry = 0; iTry < MAXTRY; iTry++) {
3473
3474	upperTheta = 1.0e50;
3475	numberPossiblySwapped = numberColumns_;
3476	numberRemaining = 0;
3477
3478	increaseInThis = 0.0; //objective increase in this loop
3479
3480	thruThis = 0.0;
3481
3482	spare = array[iFlip];
3483	index = indices[iFlip];
3484	spare2 = array[1-iFlip];
3485	index2 = indices[1-iFlip];
3486	for (i = 0; i < interesting[iFlip]; i++) {
3487	int iSequence = index[i];
3488	double alpha = spare[i];
3489	double oldValue = dj_[iSequence];
3490	double value = oldValue - upperTheta * alpha;
3491
3492	if (alpha < 0.0) {
3493	//at upper bound
3494	if (value > newTolerance) {
3495	if (-alpha >= acceptablePivot) {
3496	upperTheta = (oldValue - newTolerance) / alpha;
3497	}
3498	}
3499	} else {
3500	// at lower bound
3501	if (value < -newTolerance) {
3502	if (alpha >= acceptablePivot) {
3503	upperTheta = (oldValue + newTolerance) / alpha;
3504	}
3505	}
3506	}
3507	}
3508	bestPivot = acceptablePivot;
3509	sequenceIn_ = -1;
3510	#ifdef DUBIOUS_WEIGHTS
3511	double bestWeight = COIN_DBL_MAX;
3512	#endif
3513	double largestPivot = acceptablePivot;
3514	// now choose largest and sum all ones which will go through
3515	//printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
3516	// Sum of bad small pivots
3517	#ifdef MORE_CAREFUL
3518	double sumBadPivots = 0.0;
3519	badSumPivots = false;
3520	#endif
3521	// Make sure upperTheta will work (-O2 and above gives problems)
3522	upperTheta *= 1.0000000001;
3523	for (i = 0; i < interesting[iFlip]; i++) {
3524	int iSequence = index[i];
3525	double alpha = spare[i];
3526	double value = dj_[iSequence] - upperTheta * alpha;
3527	double badDj = 0.0;
3528
3529	bool addToSwapped = false;
3530
3531	if (alpha < 0.0) {
3532	//at upper bound
3533	if (value >= 0.0) {
3534	addToSwapped = true;
3535	#if TRYBIAS==1
3536	badDj = -CoinMax(dj_[iSequence], 0.0);
3537	#elif TRYBIAS==2
3538	badDj = -dj_[iSequence];
3539	#else
3540	badDj = -dj_[iSequence] - dualTolerance_;
3541	#endif
3542	}
3543	} else {
3544	// at lower bound
3545	if (value <= 0.0) {
3546	addToSwapped = true;
3547	#if TRYBIAS==1
3548	badDj = CoinMin(dj_[iSequence], 0.0);
3549	#elif TRYBIAS==2
3550	badDj = dj_[iSequence];
3551	#else
3552	badDj = dj_[iSequence] - dualTolerance_;
3553	#endif
3554	}
3555	}
3556	if (!addToSwapped) {
3557	// add to list of remaining
3558	spare2[numberRemaining] = alpha;
3559	index2[numberRemaining++] = iSequence;
3560	} else {
3561	// add to list of swapped
3562	spare2[--numberPossiblySwapped] = alpha;
3563	index2[numberPossiblySwapped] = iSequence;
3564	// select if largest pivot
3565	bool take = false;
3566	double absAlpha = fabs(alpha);
3567	#ifdef DUBIOUS_WEIGHTS
3568	// User could do anything to break ties here
3569	double weight;
3570	if (dubiousWeights)
3571	weight = dubiousWeights[iSequence];
3572	else
3573	weight = 1.0;
3574	weight += randomNumberGenerator_.randomDouble() * 1.0e-2;
3575	if (absAlpha > 2.0 * bestPivot) {
3576	take = true;
3577	} else if (absAlpha > largestPivot) {
3578	// could multiply absAlpha and weight
3579	if (weight * bestPivot < bestWeight * absAlpha)
3580	take = true;
3581	}
3582	#else
3583	if (absAlpha > bestPivot)
3584	take = true;
3585	#endif
3586	#ifdef MORE_CAREFUL
3587	if (absAlpha < acceptablePivot && upperTheta < 1.0e20) {
3588	if (alpha < 0.0) {
3589	//at upper bound
3590	if (value > dualTolerance_) {
3591	double gap = upper_[iSequence] - lower_[iSequence];
3592	if (gap < 1.0e20)
3593	sumBadPivots += value * gap;
3594	else
3595	sumBadPivots += 1.0e20;
3596	//printf("bad %d alpha %g dj at upper %g\n",
3597	// iSequence,alpha,value);
3598	}
3599	} else {
3600	//at lower bound
3601	if (value < -dualTolerance_) {
3602	double gap = upper_[iSequence] - lower_[iSequence];
3603	if (gap < 1.0e20)
3604	sumBadPivots -= value * gap;
3605	else
3606	sumBadPivots += 1.0e20;
3607	//printf("bad %d alpha %g dj at lower %g\n",
3608	// iSequence,alpha,value);
3609	}
3610	}
3611	}
3612	#endif
3613	#ifdef FORCE_FOLLOW
3614	if (iSequence == force_in) {
3615	printf("taking %d - alpha %g best %g\n", force_in, absAlpha, largestPivot);
3616	take = true;
3617	}
3618	#endif
3619	if (take) {
3620	sequenceIn_ = numberPossiblySwapped;
3621	bestPivot = absAlpha;
3622	theta_ = dj_[iSequence] / alpha;
3623	largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
3624	#ifdef DUBIOUS_WEIGHTS
3625	bestWeight = weight;
3626	#endif
3627	//printf(" taken seq %d alpha %g weight %d\n",
3628	// iSequence,absAlpha,dubiousWeights[iSequence]);
3629	} else {
3630	//printf(" not taken seq %d alpha %g weight %d\n",
3631	// iSequence,absAlpha,dubiousWeights[iSequence]);
3632	}
3633	double range = upper_[iSequence] - lower_[iSequence];
3634	thruThis += range * fabs(alpha);
3635	increaseInThis += badDj * range;
3636	}
3637	}
3638	marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3639	marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3640	#ifdef MORE_CAREFUL
3641	// If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3642	if (sumBadPivots > 1.0e4) {
3643	if (handler_->logLevel() > 1)
3644	printf("maybe forcing re-factorization - sum %g %d pivots\n", sumBadPivots,
3645	factorization_->pivots());
3646	if(factorization_->pivots() > 3) {
3647	badSumPivots = true;
3648	break;
3649	}
3650	}
3651	#endif
3652	swapped[1-iFlip] = numberPossiblySwapped;
3653	interesting[1-iFlip] = numberRemaining;
3654	// If we stop now this will be increase in objective (I think)
3655	double increase = (fabs(dualOut_) - totalThru) * theta_;
3656	increase += increaseInObjective;
3657	if (theta_ < 0.0)
3658	thruThis += fabs(dualOut_); // force using this one
3659	if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
3660	// back
3661	// We may need to be more careful - we could do by
3662	// switch so we always do fine grained?
3663	bestPivot = 0.0;
3664	} else {
3665	// add in
3666	totalThru += thruThis;
3667	increaseInObjective += increaseInThis;
3668	}
3669	if (bestPivot < 0.1 * bestEverPivot &&
3670	bestEverPivot > 1.0e-6 &&
3671	(bestPivot < 1.0e-3 \|\| totalThru * 2.0 > fabs(dualOut_))) {
3672	// back to previous one
3673	sequenceIn_ = lastSequence;
3674	// swap regions
3675	iFlip = 1 - iFlip;
3676	break;
3677	} else if (sequenceIn_ == -1 && upperTheta > largeValue_) {
3678	if (lastPivot > acceptablePivot) {
3679	// back to previous one
3680	sequenceIn_ = lastSequence;
3681	// swap regions
3682	iFlip = 1 - iFlip;
3683	} else {
3684	// can only get here if all pivots too small
3685	}
3686	break;
3687	} else if (totalThru >= fabs(dualOut_)) {
3688	modifyCosts = true; // fine grain - we can modify costs
3689	break; // no point trying another loop
3690	} else {
3691	lastSequence = sequenceIn_;
3692	if (bestPivot > bestEverPivot)
3693	bestEverPivot = bestPivot;
3694	iFlip = 1 - iFlip;
3695	modifyCosts = true; // fine grain - we can modify costs
3696	}
3697	}
3698	if (iTry == MAXTRY)
3699	iFlip = 1 - iFlip; // flip back
3700	break;
3701	} else {
3702	// skip this lot
3703	if (bestPivot > 1.0e-3 \|\| bestPivot > bestEverPivot) {
3704	bestEverPivot = bestPivot;
3705	lastSequence = bestSequence;
3706	} else {
3707	// keep old swapped
3708	CoinMemcpyN(array[iFlip] + swapped[iFlip],
3709	numberColumns_ - swapped[iFlip], array[1-iFlip] + swapped[iFlip]);
3710	CoinMemcpyN(indices[iFlip] + swapped[iFlip],
3711	numberColumns_ - swapped[iFlip], indices[1-iFlip] + swapped[iFlip]);
3712	marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], swapped[iFlip]);
3713	swapped[1-iFlip] = swapped[iFlip];
3714	}
3715	increaseInObjective += increaseInThis;
3716	iFlip = 1 - iFlip; // swap regions
3717	tentativeTheta = 2.0 * upperTheta;
3718	totalThru += thruThis;
3719	}
3720	}
3721
3722	// can get here without sequenceIn_ set but with lastSequence
3723	if (sequenceIn_ < 0 && lastSequence >= 0) {
3724	// back to previous one
3725	sequenceIn_ = lastSequence;
3726	// swap regions
3727	iFlip = 1 - iFlip;
3728	}
3729
3730	#define MINIMUMTHETA 1.0e-18
3731	// Movement should be minimum for anti-degeneracy - unless
3732	// fixed variable out
3733	double minimumTheta;
3734	if (upperOut_ > lowerOut_)
3735	minimumTheta = MINIMUMTHETA;
3736	else
3737	minimumTheta = 0.0;
3738	if (sequenceIn_ >= 0) {
3739	// at this stage sequenceIn_ is just pointer into index array
3740	// flip just so we can use iFlip
3741	iFlip = 1 - iFlip;
3742	spare = array[iFlip];
3743	index = indices[iFlip];
3744	double oldValue;
3745	alpha_ = spare[sequenceIn_];
3746	sequenceIn_ = indices[iFlip][sequenceIn_];
3747	oldValue = dj_[sequenceIn_];
3748	theta_ = CoinMax(oldValue / alpha_, 0.0);
3749	if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
3750	// can't pivot to zero
3751	#if 0
3752	if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) {
3753	theta_ = minimumTheta;
3754	} else if (oldValue - minimumTheta*alpha_ >= -newTolerance) {
3755	theta_ = minimumTheta;
3756	thisIncrease = true;
3757	} else {
3758	theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0);
3759	thisIncrease = true;
3760	}
3761	#else
3762	theta_ = minimumTheta;
3763	#endif
3764	}
3765	// may need to adjust costs so all dual feasible AND pivoted is exactly 0
3766	//int costOffset = numberRows_+numberColumns_;
3767	if (modifyCosts && !badSumPivots) {
3768	int i;
3769	for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) {
3770	int iSequence = index[i];
3771	double alpha = spare[i];
3772	double value = dj_[iSequence] - theta_ * alpha;
3773
3774	// can't be free here
3775
3776	if (alpha < 0.0) {
3777	//at upper bound
3778	if (value > dualTolerance_) {
3779	thisIncrease = true;
3780	#define MODIFYCOST 2
3781	#if MODIFYCOST
3782	// modify cost to hit new tolerance
3783	double modification = alpha * theta_ - dj_[iSequence]
3784	+ newTolerance;
3785	if ((specialOptions_&(2048 + 4096 + 16384)) != 0) {
3786	if ((specialOptions_ & 16384) != 0) {
3787	if (fabs(modification) < 1.0e-8)
3788	modification = 0.0;
3789	} else if ((specialOptions_ & 2048) != 0) {
3790	if (fabs(modification) < 1.0e-10)
3791	modification = 0.0;
3792	} else {
3793	if (fabs(modification) < 1.0e-12)
3794	modification = 0.0;
3795	}
3796	}
3797	dj_[iSequence] += modification;
3798	cost_[iSequence] += modification;
3799	if (modification)
3800	numberChanged_ ++; // Say changed costs
3801	//cost_[iSequence+costOffset] += modification; // save change
3802	#endif
3803	}
3804	} else {
3805	// at lower bound
3806	if (-value > dualTolerance_) {
3807	thisIncrease = true;
3808	#if MODIFYCOST
3809	// modify cost to hit new tolerance
3810	double modification = alpha * theta_ - dj_[iSequence]
3811	- newTolerance;
3812	//modification = CoinMax(modification,-dualTolerance_);
3813	//assert (fabs(modification)<1.0e-7);
3814	if ((specialOptions_&(2048 + 4096)) != 0) {
3815	if ((specialOptions_ & 2048) != 0) {
3816	if (fabs(modification) < 1.0e-10)
3817	modification = 0.0;
3818	} else {
3819	if (fabs(modification) < 1.0e-12)
3820	modification = 0.0;
3821	}
3822	}
3823	dj_[iSequence] += modification;
3824	cost_[iSequence] += modification;
3825	if (modification)
3826	numberChanged_ ++; // Say changed costs
3827	//cost_[iSequence+costOffset] += modification; // save change
3828	#endif
3829	}
3830	}
3831	}
3832	}
3833	}
3834	}
3835
3836	#ifdef MORE_CAREFUL
3837	// If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3838	if ((badSumPivots \|\|
3839	fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) {
3840	if (handler_->logLevel() > 1)
3841	printf("forcing re-factorization\n");
3842	sequenceIn_ = -1;
3843	}
3844	#endif
3845	if (sequenceIn_ >= 0) {
3846	lowerIn_ = lower_[sequenceIn_];
3847	upperIn_ = upper_[sequenceIn_];
3848	valueIn_ = solution_[sequenceIn_];
3849	dualIn_ = dj_[sequenceIn_];
3850
3851	if (numberTimesOptimal_) {
3852	// can we adjust cost back closer to original
3853	//*** add coding
3854	}
3855	#if MODIFYCOST>1
3856	// modify cost to hit zero exactly
3857	// so (dualIn_+modification)==theta_*alpha_
3858	double modification = theta_ * alpha_ - dualIn_;
3859	// But should not move objective too much ??
3860	#define DONT_MOVE_OBJECTIVE
3861	#ifdef DONT_MOVE_OBJECTIVE
3862	double moveObjective = fabs(modification * solution_[sequenceIn_]);
3863	double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3);
3864	if (moveObjective > smallMove) {
3865	if (handler_->logLevel() > 1)
3866	printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
3867	modification, solution_[sequenceIn_]);
3868	modification *= smallMove / moveObjective;
3869	}
3870	#endif
3871	if (badSumPivots)
3872	modification = 0.0;
3873	if ((specialOptions_&(2048 + 4096)) != 0) {
3874	if ((specialOptions_ & 16384) != 0) {
3875	// in fast dual
3876	if (fabs(modification) < 1.0e-7)
3877	modification = 0.0;
3878	} else if ((specialOptions_ & 2048) != 0) {
3879	if (fabs(modification) < 1.0e-10)
3880	modification = 0.0;
3881	} else {
3882	if (fabs(modification) < 1.0e-12)
3883	modification = 0.0;
3884	}
3885	}
3886	dualIn_ += modification;
3887	dj_[sequenceIn_] = dualIn_;
3888	cost_[sequenceIn_] += modification;
3889	if (modification)
3890	numberChanged_ ++; // Say changed costs
3891	//int costOffset = numberRows_+numberColumns_;
3892	//cost_[sequenceIn_+costOffset] += modification; // save change
3893	//assert (fabs(modification)<1.0e-6);
3894	#ifdef CLP_DEBUG
3895	if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
3896	printf("exact %d new cost %g, change %g\n", sequenceIn_,
3897	cost_[sequenceIn_], modification);
3898	#endif
3899	#endif
3900
3901	if (alpha_ < 0.0) {
3902	// as if from upper bound
3903	directionIn_ = -1;
3904	upperIn_ = valueIn_;
3905	} else {
3906	// as if from lower bound
3907	directionIn_ = 1;
3908	lowerIn_ = valueIn_;
3909	}
3910	} else {
3911	// no pivot
3912	bestPossible = 0.0;
3913	alpha_ = 0.0;
3914	}
3915	//if (thisIncrease)
3916	//dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3917
3918	// clear arrays
3919
3920	for (i = 0; i < 2; i++) {
3921	CoinZeroN(array[i], marker[i][0]);
3922	CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]);
3923	}
3924	return bestPossible;
3925	}
3926	#ifdef CLP_ALL_ONE_FILE
3927	#undef MAXTRY
3928	#endif
3929	/* Checks if tentative optimal actually means unbounded
3930	Returns -3 if not, 2 if is unbounded */
3931	int
3932	ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3933	CoinIndexedVector * spare,
3934	double changeCost)
3935	{
3936	int status = 2; // say unbounded
3937	factorization_->updateColumn(spare, ray);
3938	// get reduced cost
3939	int i;
3940	int number = ray->getNumElements();
3941	int * index = ray->getIndices();
3942	double * array = ray->denseVector();
3943	for (i = 0; i < number; i++) {
3944	int iRow = index[i];
3945	int iPivot = pivotVariable_[iRow];
3946	changeCost -= cost(iPivot) * array[iRow];
3947	}
3948	double way;
3949	if (changeCost > 0.0) {
3950	//try going down
3951	way = 1.0;
3952	} else if (changeCost < 0.0) {
3953	//try going up
3954	way = -1.0;
3955	} else {
3956	#ifdef CLP_DEBUG
3957	printf("can't decide on up or down\n");
3958	#endif
3959	way = 0.0;
3960	status = -3;
3961	}
3962	double movement = 1.0e10 * way; // some largish number
3963	double zeroTolerance = 1.0e-14 * dualBound_;
3964	for (i = 0; i < number; i++) {
3965	int iRow = index[i];
3966	int iPivot = pivotVariable_[iRow];
3967	double arrayValue = array[iRow];
3968	if (fabs(arrayValue) < zeroTolerance)
3969	arrayValue = 0.0;
3970	double newValue = solution(iPivot) + movement * arrayValue;
3971	if (newValue > upper(iPivot) + primalTolerance_ \|\|
3972	newValue < lower(iPivot) - primalTolerance_)
3973	status = -3; // not unbounded
3974	}
3975	if (status == 2) {
3976	// create ray
3977	delete [] ray_;
3978	ray_ = new double [numberColumns_];
3979	CoinZeroN(ray_, numberColumns_);
3980	for (i = 0; i < number; i++) {
3981	int iRow = index[i];
3982	int iPivot = pivotVariable_[iRow];
3983	double arrayValue = array[iRow];
3984	if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
3985	ray_[iPivot] = way * array[iRow];
3986	}
3987	}
3988	ray->clear();
3989	return status;
3990	}
3991	//static int count_alpha=0;
3992	/* Checks if finished. Updates status */
3993	void
3994	ClpSimplexDual::statusOfProblemInDual(int & lastCleaned, int type,
3995	double * givenDuals, ClpDataSave & saveData,
3996	int ifValuesPass)
3997	{
3998	#ifdef CLP_INVESTIGATE_SERIAL
3999	if (z_thinks > 0 && z_thinks < 2)
4000	z_thinks += 2;
4001	#endif
4002	bool arraysNotCreated = (type==0);
4003	// If lots of iterations then adjust costs if large ones
4004	if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) {
4005	double largest = 0.0;
4006	for (int i = 0; i < numberRows_; i++) {
4007	int iColumn = pivotVariable_[i];
4008	largest = CoinMax(largest, fabs(cost_[iColumn]));
4009	}
4010	if (largest > 1.0e6) {
4011	objectiveScale_ = 1.0e6 / largest;
4012	for (int i = 0; i < numberRows_ + numberColumns_; i++)
4013	cost_[i] *= objectiveScale_;
4014	}
4015	}
4016	int numberPivots = factorization_->pivots();
4017	double realDualInfeasibilities = 0.0;
4018	if (type == 2) {
4019	if (alphaAccuracy_ != -1.0)
4020	alphaAccuracy_ = -2.0;
4021	// trouble - restore solution
4022	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4023	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4024	numberRows_, rowActivityWork_);
4025	CoinMemcpyN(savedSolution_ ,
4026	numberColumns_, columnActivityWork_);
4027	// restore extra stuff
4028	int dummy;
4029	matrix_->generalExpanded(this, 6, dummy);
4030	forceFactorization_ = 1; // a bit drastic but ..
4031	changeMade_++; // say something changed
4032	// get correct bounds on all variables
4033	resetFakeBounds(0);
4034	}
4035	int tentativeStatus = problemStatus_;
4036	double changeCost;
4037	bool unflagVariables = true;
4038	bool weightsSaved = false;
4039	bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_;
4040	int dontFactorizePivots = dontFactorizePivots_;
4041	if (type == 3) {
4042	type = 1;
4043	dontFactorizePivots = 1;
4044	}
4045	if (alphaAccuracy_ < 0.0 \|\| !numberPivots \|\| alphaAccuracy_ > 1.0e4 \|\| numberPivots > 20) {
4046	if (problemStatus_ > -3 \|\| numberPivots > dontFactorizePivots) {
4047	// factorize
4048	// later on we will need to recover from singularities
4049	// also we could skip if first time
4050	// save dual weights
4051	dualRowPivot_->saveWeights(this, 1);
4052	weightsSaved = true;
4053	if (type) {
4054	// is factorization okay?
4055	if (internalFactorize(1)) {
4056	// no - restore previous basis
4057	unflagVariables = false;
4058	assert (type == 1);
4059	changeMade_++; // say something changed
4060	// Keep any flagged variables
4061	int i;
4062	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4063	if (flagged(i))
4064	saveStatus_[i] \|= 64; //say flagged
4065	}
4066	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4067	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4068	numberRows_, rowActivityWork_);
4069	CoinMemcpyN(savedSolution_ ,
4070	numberColumns_, columnActivityWork_);
4071	// restore extra stuff
4072	int dummy;
4073	matrix_->generalExpanded(this, 6, dummy);
4074	// get correct bounds on all variables
4075	resetFakeBounds(1);
4076	// need to reject something
4077	char x = isColumn(sequenceOut_) ? 'C' : 'R';
4078	handler_->message(CLP_SIMPLEX_FLAG, messages_)
4079	<< x << sequenceWithin(sequenceOut_)
4080	<< CoinMessageEol;
4081	#ifdef COIN_DEVELOP
4082	printf("flag d\n");
4083	#endif
4084	setFlagged(sequenceOut_);
4085	progress_.clearBadTimes();
4086
4087	// Go to safe
4088	factorization_->pivotTolerance(0.99);
4089	forceFactorization_ = 1; // a bit drastic but ..
4090	type = 2;
4091	//assert (internalFactorize(1)==0);
4092	if (internalFactorize(1)) {
4093	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4094	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4095	numberRows_, rowActivityWork_);
4096	CoinMemcpyN(savedSolution_ ,
4097	numberColumns_, columnActivityWork_);
4098	// restore extra stuff
4099	int dummy;
4100	matrix_->generalExpanded(this, 6, dummy);
4101	// debug
4102	int returnCode = internalFactorize(1);
4103	while (returnCode) {
4104	// ouch
4105	// switch off dense
4106	int saveDense = factorization_->denseThreshold();
4107	factorization_->setDenseThreshold(0);
4108	// Go to safe
4109	factorization_->pivotTolerance(0.99);
4110	// make sure will do safe factorization
4111	pivotVariable_[0] = -1;
4112	returnCode = internalFactorize(2);
4113	factorization_->setDenseThreshold(saveDense);
4114	}
4115	// get correct bounds on all variables
4116	resetFakeBounds(1);
4117	}
4118	}
4119	}
4120	if (problemStatus_ != -4 \|\| numberPivots > 10)
4121	problemStatus_ = -3;
4122	}
4123	} else {
4124	//printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
4125	//count_alpha++;
4126	//if ((count_alpha%5000)==0)
4127	//printf("count alpha %d\n",count_alpha);
4128	}
4129	// at this stage status is -3 or -4 if looks infeasible
4130	// get primal and dual solutions
4131	#if 0
4132	{
4133	int numberTotal = numberRows_ + numberColumns_;
4134	double * saveSol = CoinCopyOfArray(solution_, numberTotal);
4135	double * saveDj = CoinCopyOfArray(dj_, numberTotal);
4136	double tolerance = type ? 1.0e-4 : 1.0e-8;
4137	// always if values pass
4138	double saveObj = objectiveValue_;
4139	double sumPrimal = sumPrimalInfeasibilities_;
4140	int numberPrimal = numberPrimalInfeasibilities_;
4141	double sumDual = sumDualInfeasibilities_;
4142	int numberDual = numberDualInfeasibilities_;
4143	gutsOfSolution(givenDuals, NULL);
4144	int j;
4145	double largestPrimal = tolerance;
4146	int iPrimal = -1;
4147	for (j = 0; j < numberTotal; j++) {
4148	double difference = solution_[j] - saveSol[j];
4149	if (fabs(difference) > largestPrimal) {
4150	iPrimal = j;
4151	largestPrimal = fabs(difference);
4152	}
4153	}
4154	double largestDual = tolerance;
4155	int iDual = -1;
4156	for (j = 0; j < numberTotal; j++) {
4157	double difference = dj_[j] - saveDj[j];
4158	if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
4159	iDual = j;
4160	largestDual = fabs(difference);
4161	}
4162	}
4163	if (!type) {
4164	if (fabs(saveObj - objectiveValue_) > 1.0e-5 \|\|
4165	numberPrimal != numberPrimalInfeasibilities_ \|\| numberPrimal != 1 \|\|
4166	fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 \|\| iPrimal >= 0 \|\|
4167	numberDual != numberDualInfeasibilities_ \|\| numberDual != 0 \|\|
4168	fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 \|\| iDual >= 0)
4169	printf("type %d its %d pivots %d primal n(%d,%d) s(%g,%g) diff(%g,%d) dual n(%d,%d) s(%g,%g) diff(%g,%d) obj(%g,%g)\n",
4170	type, numberIterations_, numberPivots,
4171	numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4172	largestPrimal, iPrimal,
4173	numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4174	largestDual, iDual,
4175	saveObj, objectiveValue_);
4176	} else {
4177	if (fabs(saveObj - objectiveValue_) > 1.0e-5 \|\|
4178	numberPrimalInfeasibilities_ \|\| iPrimal >= 0 \|\|
4179	numberDualInfeasibilities_ \|\| iDual >= 0)
4180	printf("type %d its %d pivots %d primal n(%d,%d) s(%g,%g) diff(%g,%d) dual n(%d,%d) s(%g,%g) diff(%g,%d) obj(%g,%g)\n",
4181	type, numberIterations_, numberPivots,
4182	numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4183	largestPrimal, iPrimal,
4184	numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4185	largestDual, iDual,
4186	saveObj, objectiveValue_);
4187	}
4188	delete [] saveSol;
4189	delete [] saveDj;
4190	}
4191	#else
4192	if (type \|\| ifValuesPass)
4193	gutsOfSolution(givenDuals, NULL);
4194	#endif
4195	// If bad accuracy treat as singular
4196	if ((largestPrimalError_ > 1.0e15 \|\| largestDualError_ > 1.0e15) && numberIterations_) {
4197	// restore previous basis
4198	unflagVariables = false;
4199	changeMade_++; // say something changed
4200	// Keep any flagged variables
4201	int i;
4202	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4203	if (flagged(i))
4204	saveStatus_[i] \|= 64; //say flagged
4205	}
4206	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4207	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4208	numberRows_, rowActivityWork_);
4209	CoinMemcpyN(savedSolution_ ,
4210	numberColumns_, columnActivityWork_);
4211	// restore extra stuff
4212	int dummy;
4213	matrix_->generalExpanded(this, 6, dummy);
4214	// get correct bounds on all variables
4215	resetFakeBounds(1);
4216	// need to reject something
4217	char x = isColumn(sequenceOut_) ? 'C' : 'R';
4218	handler_->message(CLP_SIMPLEX_FLAG, messages_)
4219	<< x << sequenceWithin(sequenceOut_)
4220	<< CoinMessageEol;
4221	#ifdef COIN_DEVELOP
4222	printf("flag e\n");
4223	#endif
4224	setFlagged(sequenceOut_);
4225	progress_.clearBadTimes();
4226
4227	// Go to safer
4228	double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99);
4229	factorization_->pivotTolerance(newTolerance);
4230	forceFactorization_ = 1; // a bit drastic but ..
4231	if (alphaAccuracy_ != -1.0)
4232	alphaAccuracy_ = -2.0;
4233	type = 2;
4234	//assert (internalFactorize(1)==0);
4235	if (internalFactorize(1)) {
4236	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4237	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4238	numberRows_, rowActivityWork_);
4239	CoinMemcpyN(savedSolution_ ,
4240	numberColumns_, columnActivityWork_);
4241	// restore extra stuff
4242	int dummy;
4243	matrix_->generalExpanded(this, 6, dummy);
4244	// debug
4245	int returnCode = internalFactorize(1);
4246	while (returnCode) {
4247	// ouch
4248	// switch off dense
4249	int saveDense = factorization_->denseThreshold();
4250	factorization_->setDenseThreshold(0);
4251	// Go to safe
4252	factorization_->pivotTolerance(0.99);
4253	// make sure will do safe factorization
4254	pivotVariable_[0] = -1;
4255	returnCode = internalFactorize(2);
4256	factorization_->setDenseThreshold(saveDense);
4257	}
4258	// get correct bounds on all variables
4259	resetFakeBounds(1);
4260	}
4261	// get primal and dual solutions
4262	gutsOfSolution(givenDuals, NULL);
4263	} else if (goodAccuracy()) {
4264	// Can reduce tolerance
4265	double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
4266	factorization_->pivotTolerance(newTolerance);
4267	}
4268	bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
4269	objectiveValue_ - bestPossibleImprovement_);
4270	bool reallyBadProblems = false;
4271	// Double check infeasibility if no action
4272	if (progress_.lastIterationNumber(0) == numberIterations_) {
4273	if (dualRowPivot_->looksOptimal()) {
4274	numberPrimalInfeasibilities_ = 0;
4275	sumPrimalInfeasibilities_ = 0.0;
4276	}
4277	#if 1
4278	} else {
4279	double thisObj = objectiveValue_ - bestPossibleImprovement_;
4280	#ifdef CLP_INVESTIGATE
4281	assert (bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100);
4282	if (bestPossibleImprovement_)
4283	printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_,
4284	thisObj);
4285	#endif
4286	double lastObj = progress_.lastObjective(0);
4287	#ifndef NDEBUG
4288	#ifdef COIN_DEVELOP
4289	resetFakeBounds(-1);
4290	#endif
4291	#endif
4292	#ifdef CLP_REPORT_PROGRESS
4293	ixxxxxx++;
4294	if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4295	char temp[20];
4296	sprintf(temp, "sol%d.out", ixxxxxx);
4297	printf("sol%d.out\n", ixxxxxx);
4298	FILE * fp = fopen(temp, "w");
4299	int nTotal = numberRows_ + numberColumns_;
4300	for (int i = 0; i < nTotal; i++)
4301	fprintf(fp, "%d %d %g %g %g %g %g\n",
4302	i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4303	fclose(fp);
4304	}
4305	#endif
4306	if(!ifValuesPass && firstFree_ < 0) {
4307	double testTol = 5.0e-3;
4308	if (progress_.timesFlagged() > 10) {
4309	testTol *= pow(2.0, progress_.timesFlagged() - 8);
4310	} else if (progress_.timesFlagged() > 5) {
4311	testTol *= 5.0;
4312	}
4313	if (lastObj > thisObj +
4314	testTol*(fabs(thisObj) + fabs(lastObj)) + testTol) {
4315	int maxFactor = factorization_->maximumPivots();
4316	if ((specialOptions_ & 1048576) == 0) {
4317	if (progress_.timesFlagged() > 10)
4318	progress_.incrementReallyBadTimes();
4319	if (maxFactor > 10 - 9) {
4320	#ifdef COIN_DEVELOP
4321	printf("lastobj %g thisobj %g\n", lastObj, thisObj);
4322	#endif
4323	//if (forceFactorization_<0)
4324	//forceFactorization_= maxFactor;
4325	//forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4326	if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 &&
4327	largestPrimalError_ < 1.0e2) {
4328	// Just save costs
4329	// save extra copy of cost_
4330	int nTotal = numberRows_ + numberColumns_;
4331	double * temp = new double [2*nTotal];
4332	memcpy(temp, cost_, nTotal * sizeof(double));
4333	memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4334	delete [] cost_;
4335	cost_ = temp;
4336	objectiveWork_ = cost_;
4337	rowObjectiveWork_ = cost_ + numberColumns_;
4338	progressFlag_ \|= 4;
4339	} else {
4340	forceFactorization_ = 1;
4341	#ifdef COIN_DEVELOP
4342	printf("Reducing factorization frequency - bad backwards\n");
4343	#endif
4344	#if 1
4345	unflagVariables = false;
4346	changeMade_++; // say something changed
4347	int nTotal = numberRows_ + numberColumns_;
4348	CoinMemcpyN(saveStatus_, nTotal, status_);
4349	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4350	numberRows_, rowActivityWork_);
4351	CoinMemcpyN(savedSolution_ ,
4352	numberColumns_, columnActivityWork_);
4353	if ((progressFlag_ & 4) == 0) {
4354	// save extra copy of cost_
4355	double * temp = new double [2*nTotal];
4356	memcpy(temp, cost_, nTotal * sizeof(double));
4357	memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4358	delete [] cost_;
4359	cost_ = temp;
4360	objectiveWork_ = cost_;
4361	rowObjectiveWork_ = cost_ + numberColumns_;
4362	progressFlag_ \|= 4;
4363	} else {
4364	memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double));
4365	}
4366	// restore extra stuff
4367	int dummy;
4368	matrix_->generalExpanded(this, 6, dummy);
4369	double pivotTolerance = factorization_->pivotTolerance();
4370	if(pivotTolerance < 0.2)
4371	factorization_->pivotTolerance(0.2);
4372	else if(progress_.timesFlagged() > 2)
4373	factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99));
4374	if (alphaAccuracy_ != -1.0)
4375	alphaAccuracy_ = -2.0;
4376	if (internalFactorize(1)) {
4377	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4378	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4379	numberRows_, rowActivityWork_);
4380	CoinMemcpyN(savedSolution_ ,
4381	numberColumns_, columnActivityWork_);
4382	// restore extra stuff
4383	int dummy;
4384	matrix_->generalExpanded(this, 6, dummy);
4385	// debug
4386	int returnCode = internalFactorize(1);
4387	while (returnCode) {
4388	// ouch
4389	// switch off dense
4390	int saveDense = factorization_->denseThreshold();
4391	factorization_->setDenseThreshold(0);
4392	// Go to safe
4393	factorization_->pivotTolerance(0.99);
4394	// make sure will do safe factorization
4395	pivotVariable_[0] = -1;
4396	returnCode = internalFactorize(2);
4397	factorization_->setDenseThreshold(saveDense);
4398	}
4399	}
4400	resetFakeBounds(0);
4401	type = 2; // so will restore weights
4402	// get primal and dual solutions
4403	gutsOfSolution(givenDuals, NULL);
4404	if (numberPivots < 2) {
4405	// need to reject something
4406	char x = isColumn(sequenceOut_) ? 'C' : 'R';
4407	handler_->message(CLP_SIMPLEX_FLAG, messages_)
4408	<< x << sequenceWithin(sequenceOut_)
4409	<< CoinMessageEol;
4410	#ifdef COIN_DEVELOP
4411	printf("flag d\n");
4412	#endif
4413	setFlagged(sequenceOut_);
4414	progress_.clearBadTimes();
4415	progress_.incrementTimesFlagged();
4416	}
4417	if (numberPivots < 10)
4418	reallyBadProblems = true;
4419	#ifdef COIN_DEVELOP
4420	printf("obj now %g\n", objectiveValue_);
4421	#endif
4422	progress_.modifyObjective(objectiveValue_
4423	- bestPossibleImprovement_);
4424	#endif
4425	}
4426	}
4427	} else {
4428	// in fast dual give up
4429	#ifdef COIN_DEVELOP
4430	printf("In fast dual?\n");
4431	#endif
4432	problemStatus_ = 3;
4433	}
4434	} else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) {
4435	numberTimesOptimal_ = 0;
4436	}
4437	}
4438	#endif
4439	}
4440	// Up tolerance if looks a bit odd
4441	if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
4442	if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
4443	int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1);
4444	if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
4445	if (factorization_->pivotTolerance() < 0.9) {
4446	// up tolerance
4447	factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
4448	//printf("tol now %g\n",factorization_->pivotTolerance());
4449	progress_.clearIterationNumbers();
4450	}
4451	}
4452	}
4453	}
4454	// Check if looping
4455	int loop;
4456	if (!givenDuals && type != 2)
4457	loop = progress_.looping();
4458	else
4459	loop = -1;
4460	if (progress_.reallyBadTimes() > 10) {
4461	problemStatus_ = 10; // instead - try other algorithm
4462	#if COIN_DEVELOP>2
4463	printf("returning at %d\n", __LINE__);
4464	#endif
4465	}
4466	int situationChanged = 0;
4467	if (loop >= 0) {
4468	problemStatus_ = loop; //exit if in loop
4469	if (!problemStatus_) {
4470	// declaring victory
4471	numberPrimalInfeasibilities_ = 0;
4472	sumPrimalInfeasibilities_ = 0.0;
4473	} else {
4474	problemStatus_ = 10; // instead - try other algorithm
4475	#if COIN_DEVELOP>2
4476	printf("returning at %d\n", __LINE__);
4477	#endif
4478	}
4479	return;
4480	} else if (loop < -1) {
4481	// something may have changed
4482	gutsOfSolution(NULL, NULL);
4483	situationChanged = 1;
4484	}
4485	// really for free variables in
4486	if((progressFlag_ & 2) != 0) {
4487	situationChanged = 2;
4488	}
4489	progressFlag_ &= (~3); //reset progress flag
4490	if ((progressFlag_ & 4) != 0) {
4491	// save copy of cost_
4492	int nTotal = numberRows_ + numberColumns_;
4493	memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double));
4494	}
4495	/*if (!numberIterations_&&sumDualInfeasibilities_)
4496	printf("OBJ %g sumPinf %g sumDinf %g\n",
4497	objectiveValue(),sumPrimalInfeasibilities_,
4498	sumDualInfeasibilities_);*/
4499	// mark as having gone optimal if looks like it
4500	if (!numberPrimalInfeasibilities_&&
4501	!numberDualInfeasibilities_)
4502	progressFlag_ \|= 8;
4503	if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
4504	handler_->message(CLP_SIMPLEX_STATUS, messages_)
4505	<< numberIterations_ << objectiveValue();
4506	handler_->printing(sumPrimalInfeasibilities_ > 0.0)
4507	<< sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
4508	handler_->printing(sumDualInfeasibilities_ > 0.0)
4509	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
4510	handler_->printing(numberDualInfeasibilitiesWithoutFree_
4511	< numberDualInfeasibilities_)
4512	<< numberDualInfeasibilitiesWithoutFree_;
4513	handler_->message() << CoinMessageEol;
4514	}
4515	#if 0
4516	printf("IT %d %g %g(%d) %g(%d)\n",
4517	numberIterations_, objectiveValue(),
4518	sumPrimalInfeasibilities_, numberPrimalInfeasibilities_,
4519	sumDualInfeasibilities_, numberDualInfeasibilities_);
4520	#endif
4521	double approximateObjective = objectiveValue_;
4522	#ifdef CLP_REPORT_PROGRESS
4523	if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4524	char temp[20];
4525	sprintf(temp, "x_sol%d.out", ixxxxxx);
4526	FILE * fp = fopen(temp, "w");
4527	int nTotal = numberRows_ + numberColumns_;
4528	for (int i = 0; i < nTotal; i++)
4529	fprintf(fp, "%d %d %g %g %g %g %g\n",
4530	i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4531	fclose(fp);
4532	if (ixxxxxx == ixxyyyy)
4533	exit(6);
4534	}
4535	#endif
4536	realDualInfeasibilities = sumDualInfeasibilities_;
4537	double saveTolerance = dualTolerance_;
4538	// If we need to carry on cleaning variables
4539	if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
4540	for (int iRow = 0; iRow < numberRows_; iRow++) {
4541	int iPivot = pivotVariable_[iRow];
4542	if (!flagged(iPivot) && pivoted(iPivot)) {
4543	// carry on
4544	numberPrimalInfeasibilities_ = -1;
4545	sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4546	sumPrimalInfeasibilities_ = 1.0;
4547	break;
4548	}
4549	}
4550	}
4551	/* If we are primal feasible and any dual infeasibilities are on
4552	free variables then it is better to go to primal */
4553	if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
4554	numberDualInfeasibilities_)
4555	problemStatus_ = 10;
4556	// dual bound coming in
4557	//double saveDualBound = dualBound_;
4558	bool needCleanFake = false;
4559	while (problemStatus_ <= -3) {
4560	int cleanDuals = 0;
4561	if (situationChanged != 0)
4562	cleanDuals = 1;
4563	int numberChangedBounds = 0;
4564	int doOriginalTolerance = 0;
4565	if ( lastCleaned == numberIterations_)
4566	doOriginalTolerance = 1;
4567	// check optimal
4568	// give code benefit of doubt
4569	if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
4570	sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4571	// say optimal (with these bounds etc)
4572	numberDualInfeasibilities_ = 0;
4573	sumDualInfeasibilities_ = 0.0;
4574	numberPrimalInfeasibilities_ = 0;
4575	sumPrimalInfeasibilities_ = 0.0;
4576	}
4577	//if (dualFeasible()\|\|problemStatus_==-4\|\|(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4578	if (dualFeasible() \|\| problemStatus_ == -4) {
4579	progress_.modifyObjective(objectiveValue_
4580	- bestPossibleImprovement_);
4581	#ifdef COIN_DEVELOP
4582	if (sumDualInfeasibilities_ \|\| bestPossibleImprovement_)
4583	printf("improve %g dualinf %g -> %g\n",
4584	bestPossibleImprovement_, sumDualInfeasibilities_,
4585	sumDualInfeasibilities_ * dualBound_);
4586	#endif
4587	// see if cutoff reached
4588	double limit = 0.0;
4589	getDblParam(ClpDualObjectiveLimit, limit);
4590	#if 0
4591	if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
4592	limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) {
4593	//looks infeasible on objective
4594	if (perturbation_ == 101) {
4595	cleanDuals = 1;
4596	// Save costs
4597	int numberTotal = numberRows_ + numberColumns_;
4598	double * saveCost = CoinCopyOfArray(cost_, numberTotal);
4599	// make sure fake bounds are back
4600	changeBounds(1, NULL, changeCost);
4601	createRim4(false);
4602	// make sure duals are current
4603	computeDuals(givenDuals);
4604	checkDualSolution();
4605	if(objectiveValue()*optimizationDirection_ >
4606	limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) {
4607	perturbation_ = 102; // stop any perturbations
4608	printf("cutoff test succeeded\n");
4609	} else {
4610	printf("cutoff test failed\n");
4611	// put back
4612	memcpy(cost_, saveCost, numberTotal * sizeof(double));
4613	// make sure duals are current
4614	computeDuals(givenDuals);
4615	checkDualSolution();
4616	progress_.modifyObjective(-COIN_DBL_MAX);
4617	problemStatus_ = -1;
4618	}
4619	delete [] saveCost;
4620	}
4621	}
4622	#endif
4623	if (primalFeasible() && !givenDuals) {
4624	// may be optimal - or may be bounds are wrong
4625	handler_->message(CLP_DUAL_BOUNDS, messages_)
4626	<< dualBound_
4627	<< CoinMessageEol;
4628	// save solution in case unbounded
4629	double * saveColumnSolution = NULL;
4630	double * saveRowSolution = NULL;
4631	bool inCbc = (specialOptions_ & (0x01000000 \| 16384)) != 0;
4632	if (!inCbc) {
4633	saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_);
4634	saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_);
4635	}
4636	numberChangedBounds = changeBounds(0, rowArray_[3], changeCost);
4637	if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
4638	//looks optimal - do we need to reset tolerance
4639	if (perturbation_ == 101) {
4640	perturbation_ = 102; // stop any perturbations
4641	cleanDuals = 1;
4642	// make sure fake bounds are back
4643	//computeObjectiveValue();
4644	changeBounds(1, NULL, changeCost);
4645	//computeObjectiveValue();
4646	createRim4(false);
4647	// make sure duals are current
4648	computeDuals(givenDuals);
4649	checkDualSolution();
4650	progress_.modifyObjective(-COIN_DBL_MAX);
4651	#define DUAL_TRY_FASTER
4652	#ifdef DUAL_TRY_FASTER
4653	if (numberDualInfeasibilities_) {
4654	#endif
4655	numberChanged_ = 1; // force something to happen
4656	lastCleaned = numberIterations_ - 1;
4657	#ifdef DUAL_TRY_FASTER
4658	} else {
4659	//double value = objectiveValue_;
4660	computeObjectiveValue(true);
4661	//printf("old %g new %g\n",value,objectiveValue_);
4662	//numberChanged_=1;
4663	}
4664	#endif
4665	}
4666	if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 &&
4667	(numberChanged_ \|\| (specialOptions_ & 4096) == 0)) {
4668	if ((specialOptions_&2097152)==0) {
4669	doOriginalTolerance = 2;
4670	numberTimesOptimal_++;
4671	changeMade_++; // say something changed
4672	if (numberTimesOptimal_ == 1) {
4673	dualTolerance_ = dblParam_[ClpDualTolerance];
4674	} else {
4675	if (numberTimesOptimal_ == 2) {
4676	// better to have small tolerance even if slower
4677	factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
4678	}
4679	dualTolerance_ = dblParam_[ClpDualTolerance];
4680	dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1);
4681	}
4682	cleanDuals = 2; // If nothing changed optimal else primal
4683	} else {
4684	// no cost - skip checks
4685	problemStatus_=0;
4686	}
4687	} else {
4688	problemStatus_ = 0; // optimal
4689	if (lastCleaned < numberIterations_ && numberChanged_) {
4690	handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
4691	<< CoinMessageEol;
4692	}
4693	}
4694	} else {
4695	cleanDuals = 1;
4696	if (doOriginalTolerance == 1) {
4697	// check unbounded
4698	// find a variable with bad dj
4699	int iSequence;
4700	int iChosen = -1;
4701	if (!inCbc) {
4702	double largest = 100.0 * primalTolerance_;
4703	for (iSequence = 0; iSequence < numberRows_ + numberColumns_;
4704	iSequence++) {
4705	double djValue = dj_[iSequence];
4706	double originalLo = originalLower(iSequence);
4707	double originalUp = originalUpper(iSequence);
4708	if (fabs(djValue) > fabs(largest)) {
4709	if (getStatus(iSequence) != basic) {
4710	if (djValue > 0 && originalLo < -1.0e20) {
4711	if (djValue > fabs(largest)) {
4712	largest = djValue;
4713	iChosen = iSequence;
4714	}
4715	} else if (djValue < 0 && originalUp > 1.0e20) {
4716	if (-djValue > fabs(largest)) {
4717	largest = djValue;
4718	iChosen = iSequence;
4719	}
4720	}
4721	}
4722	}
4723	}
4724	}
4725	if (iChosen >= 0) {
4726	int iSave = sequenceIn_;
4727	sequenceIn_ = iChosen;
4728	unpack(rowArray_[1]);
4729	sequenceIn_ = iSave;
4730	// if dual infeasibilities then must be free vector so add in dual
4731	if (numberDualInfeasibilities_) {
4732	if (fabs(changeCost) > 1.0e-5)
4733	COIN_DETAIL_PRINT(printf("Odd free/unbounded combo\n"));
4734	changeCost += cost_[iChosen];
4735	}
4736	problemStatus_ = checkUnbounded(rowArray_[1], rowArray_[0],
4737	changeCost);
4738	rowArray_[1]->clear();
4739	} else {
4740	problemStatus_ = -3;
4741	}
4742	if (problemStatus_ == 2 && perturbation_ == 101) {
4743	perturbation_ = 102; // stop any perturbations
4744	cleanDuals = 1;
4745	createRim4(false);
4746	progress_.modifyObjective(-COIN_DBL_MAX);
4747	problemStatus_ = -1;
4748	}
4749	if (problemStatus_ == 2) {
4750	// it is unbounded - restore solution
4751	// but first add in changes to non-basic
4752	int iColumn;
4753	double * original = columnArray_[0]->denseVector();
4754	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4755	if(getColumnStatus(iColumn) != basic)
4756	ray_[iColumn] +=
4757	saveColumnSolution[iColumn] - original[iColumn];
4758	columnActivityWork_[iColumn] = original[iColumn];
4759	}
4760	CoinMemcpyN(saveRowSolution, numberRows_,
4761	rowActivityWork_);
4762	}
4763	} else {
4764	doOriginalTolerance = 2;
4765	rowArray_[0]->clear();
4766	}
4767	}
4768	delete [] saveColumnSolution;
4769	delete [] saveRowSolution;
4770	}
4771	if (problemStatus_ == -4 \|\| problemStatus_ == -5) {
4772	// may be infeasible - or may be bounds are wrong
4773	numberChangedBounds = changeBounds(0, NULL, changeCost);
4774	needCleanFake = true;
4775	/* Should this be here as makes no difference to being feasible.
4776	But seems to make a difference to run times. */
4777	if (perturbation_ == 101 && 0) {
4778	perturbation_ = 102; // stop any perturbations
4779	cleanDuals = 1;
4780	numberChangedBounds = 1;
4781	// make sure fake bounds are back
4782	changeBounds(1, NULL, changeCost);
4783	needCleanFake = true;
4784	createRim4(false);
4785	progress_.modifyObjective(-COIN_DBL_MAX);
4786	}
4787	if ((numberChangedBounds <= 0 \|\| dualBound_ > 1.0e20 \|\|
4788	(largestPrimalError_ > 1.0 && dualBound_ > 1.0e17)) &&
4789	(numberPivots < 4 \|\| sumPrimalInfeasibilities_ > 1.0e-6)) {
4790	problemStatus_ = 1; // infeasible
4791	if (perturbation_ == 101) {
4792	perturbation_ = 102; // stop any perturbations
4793	//cleanDuals=1;
4794	//numberChangedBounds=1;
4795	//createRim4(false);
4796	}
4797	} else {
4798	problemStatus_ = -1; //iterate
4799	cleanDuals = 1;
4800	if (numberChangedBounds <= 0)
4801	doOriginalTolerance = 2;
4802	// and delete ray which has been created
4803	delete [] ray_;
4804	ray_ = NULL;
4805	}
4806
4807	}
4808	} else {
4809	cleanDuals = 1;
4810	}
4811	if (problemStatus_ < 0) {
4812	if (doOriginalTolerance == 2) {
4813	// put back original tolerance
4814	lastCleaned = numberIterations_;
4815	numberChanged_ = 0; // Number of variables with changed costs
4816	handler_->message(CLP_DUAL_ORIGINAL, messages_)
4817	<< CoinMessageEol;
4818	perturbation_ = 102; // stop any perturbations
4819	#if 0
4820	double * xcost = new double[numberRows_+numberColumns_];
4821	double * xlower = new double[numberRows_+numberColumns_];
4822	double * xupper = new double[numberRows_+numberColumns_];
4823	double * xdj = new double[numberRows_+numberColumns_];
4824	double * xsolution = new double[numberRows_+numberColumns_];
4825	CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4826	CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4827	CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4828	CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4829	CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4830	#endif
4831	createRim4(false);
4832	progress_.modifyObjective(-COIN_DBL_MAX);
4833	// make sure duals are current
4834	computeDuals(givenDuals);
4835	checkDualSolution();
4836	#if 0
4837	int i;
4838	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4839	if (cost_[i] != xcost[i])
4840	printf("** %d old cost %g new %g sol %g\n",
4841	i, xcost[i], cost_[i], solution_[i]);
4842	if (lower_[i] != xlower[i])
4843	printf("** %d old lower %g new %g sol %g\n",
4844	i, xlower[i], lower_[i], solution_[i]);
4845	if (upper_[i] != xupper[i])
4846	printf("** %d old upper %g new %g sol %g\n",
4847	i, xupper[i], upper_[i], solution_[i]);
4848	if (dj_[i] != xdj[i])
4849	printf("** %d old dj %g new %g sol %g\n",
4850	i, xdj[i], dj_[i], solution_[i]);
4851	if (solution_[i] != xsolution[i])
4852	printf("** %d old solution %g new %g sol %g\n",
4853	i, xsolution[i], solution_[i], solution_[i]);
4854	}
4855	//delete [] xcost;
4856	//delete [] xupper;
4857	//delete [] xlower;
4858	//delete [] xdj;
4859	//delete [] xsolution;
4860	#endif
4861	// put back bounds as they were if was optimal
4862	if (doOriginalTolerance == 2 && cleanDuals != 2) {
4863	changeMade_++; // say something changed
4864	/* We may have already changed some bounds in this function
4865	so save numberFake_ and add in.
4866
4867	Worst that can happen is that we waste a bit of time - but it must be finite.
4868	*/
4869	//int saveNumberFake = numberFake_;
4870	//resetFakeBounds(-1);
4871	changeBounds(3, NULL, changeCost);
4872	needCleanFake = true;
4873	//numberFake_ += saveNumberFake;
4874	//resetFakeBounds(-1);
4875	cleanDuals = 2;
4876	//cleanDuals=1;
4877	}
4878	#if 0
4879	//int i;
4880	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4881	if (cost_[i] != xcost[i])
4882	printf("** %d old cost %g new %g sol %g\n",
4883	i, xcost[i], cost_[i], solution_[i]);
4884	if (lower_[i] != xlower[i])
4885	printf("** %d old lower %g new %g sol %g\n",
4886	i, xlower[i], lower_[i], solution_[i]);
4887	if (upper_[i] != xupper[i])
4888	printf("** %d old upper %g new %g sol %g\n",
4889	i, xupper[i], upper_[i], solution_[i]);
4890	if (dj_[i] != xdj[i])
4891	printf("** %d old dj %g new %g sol %g\n",
4892	i, xdj[i], dj_[i], solution_[i]);
4893	if (solution_[i] != xsolution[i])
4894	printf("** %d old solution %g new %g sol %g\n",
4895	i, xsolution[i], solution_[i], solution_[i]);
4896	}
4897	delete [] xcost;
4898	delete [] xupper;
4899	delete [] xlower;
4900	delete [] xdj;
4901	delete [] xsolution;
4902	#endif
4903	}
4904	if (cleanDuals == 1 \|\| (cleanDuals == 2 && !numberDualInfeasibilities_)) {
4905	// make sure dual feasible
4906	// look at all rows and columns
4907	rowArray_[0]->clear();
4908	columnArray_[0]->clear();
4909	double objectiveChange = 0.0;
4910	double savePrimalInfeasibilities = sumPrimalInfeasibilities_;
4911	if (!numberIterations_) {
4912	int nTotal = numberRows_ + numberColumns_;
4913	if (arraysNotCreated) {
4914	// create save arrays
4915	delete [] saveStatus_;
4916	delete [] savedSolution_;
4917	saveStatus_ = new unsigned char [nTotal];
4918	savedSolution_ = new double [nTotal];
4919	arraysNotCreated = false;
4920	}
4921	// save arrays
4922	CoinMemcpyN(status_, nTotal, saveStatus_);
4923	CoinMemcpyN(rowActivityWork_,
4924	numberRows_, savedSolution_ + numberColumns_);
4925	CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
4926	}
4927	#if 0
4928	double * xcost = new double[numberRows_+numberColumns_];
4929	double * xlower = new double[numberRows_+numberColumns_];
4930	double * xupper = new double[numberRows_+numberColumns_];
4931	double * xdj = new double[numberRows_+numberColumns_];
4932	double * xsolution = new double[numberRows_+numberColumns_];
4933	CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4934	CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4935	CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4936	CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4937	CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4938	#endif
4939	if (givenDuals)
4940	dualTolerance_ = 1.0e50;
4941	updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
4942	0.0, objectiveChange, true);
4943	dualTolerance_ = saveTolerance;
4944	#if 0
4945	int i;
4946	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4947	if (cost_[i] != xcost[i])
4948	printf("** %d old cost %g new %g sol %g\n",
4949	i, xcost[i], cost_[i], solution_[i]);
4950	if (lower_[i] != xlower[i])
4951	printf("** %d old lower %g new %g sol %g\n",
4952	i, xlower[i], lower_[i], solution_[i]);
4953	if (upper_[i] != xupper[i])
4954	printf("** %d old upper %g new %g sol %g\n",
4955	i, xupper[i], upper_[i], solution_[i]);
4956	if (dj_[i] != xdj[i])
4957	printf("** %d old dj %g new %g sol %g\n",
4958	i, xdj[i], dj_[i], solution_[i]);
4959	if (solution_[i] != xsolution[i])
4960	printf("** %d old solution %g new %g sol %g\n",
4961	i, xsolution[i], solution_[i], solution_[i]);
4962	}
4963	delete [] xcost;
4964	delete [] xupper;
4965	delete [] xlower;
4966	delete [] xdj;
4967	delete [] xsolution;
4968	#endif
4969	// for now - recompute all
4970	gutsOfSolution(NULL, NULL);
4971	if (givenDuals)
4972	dualTolerance_ = 1.0e50;
4973	updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
4974	0.0, objectiveChange, true);
4975	dualTolerance_ = saveTolerance;
4976	if (!numberIterations_ && sumPrimalInfeasibilities_ >
4977	1.0e5*(savePrimalInfeasibilities+1.0e3) &&
4978	(moreSpecialOptions_ & 256) == 0) {
4979	// Use primal
4980	int nTotal = numberRows_ + numberColumns_;
4981	CoinMemcpyN(saveStatus_, nTotal, status_);
4982	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4983	numberRows_, rowActivityWork_);
4984	CoinMemcpyN(savedSolution_ ,
4985	numberColumns_, columnActivityWork_);
4986	problemStatus_ = 10;
4987	situationChanged = 0;
4988	}
4989	//assert(numberDualInfeasibilitiesWithoutFree_==0);
4990	if (numberDualInfeasibilities_) {
4991	if ((numberPrimalInfeasibilities_ \|\| numberPivots)
4992	&& problemStatus_!=10) {
4993	problemStatus_ = -1; // carry on as normal
4994	} else {
4995	problemStatus_ = 10; // try primal
4996	#if COIN_DEVELOP>1
4997	printf("returning at %d\n", __LINE__);
4998	#endif
4999	}
5000	} else if (situationChanged == 2) {
5001	problemStatus_ = -1; // carry on as normal
5002	// need to reset bounds
5003	changeBounds(3, NULL, changeCost);
5004	}
5005	situationChanged = 0;
5006	} else {
5007	// iterate
5008	if (cleanDuals != 2) {
5009	problemStatus_ = -1;
5010	} else {
5011	problemStatus_ = 10; // try primal
5012	#if COIN_DEVELOP>2
5013	printf("returning at %d\n", __LINE__);
5014	#endif
5015	}
5016	}
5017	}
5018	}
5019	// unflag all variables (we may want to wait a bit?)
5020	if ((tentativeStatus != -2 && tentativeStatus != -1) && unflagVariables) {
5021	int iRow;
5022	int numberFlagged = 0;
5023	for (iRow = 0; iRow < numberRows_; iRow++) {
5024	int iPivot = pivotVariable_[iRow];
5025	if (flagged(iPivot)) {
5026	numberFlagged++;
5027	clearFlagged(iPivot);
5028	}
5029	}
5030	#ifdef COIN_DEVELOP
5031	if (numberFlagged) {
5032	printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus,
5033	problemStatus_, numberPrimalInfeasibilities_,
5034	numberTimesOptimal_);
5035	}
5036	#endif
5037	unflagVariables = numberFlagged > 0;
5038	if (numberFlagged && !numberPivots) {
5039	/* looks like trouble as we have not done any iterations.
5040	Try changing pivot tolerance then give it a few goes and give up */
5041	if (factorization_->pivotTolerance() < 0.9) {
5042	factorization_->pivotTolerance(0.99);
5043	problemStatus_ = -1;
5044	} else if (numberTimesOptimal_ < 3) {
5045	numberTimesOptimal_++;
5046	problemStatus_ = -1;
5047	} else {
5048	unflagVariables = false;
5049	//secondaryStatus_ = 1; // and say probably infeasible
5050	if ((moreSpecialOptions_ & 256) == 0) {
5051	// try primal
5052	problemStatus_ = 10;
5053	} else {
5054	// almost certainly infeasible
5055	problemStatus_ = 1;
5056	}
5057	#if COIN_DEVELOP>1
5058	printf("returning at %d\n", __LINE__);
5059	#endif
5060	}
5061	}
5062	}
5063	if (problemStatus_ < 0) {
5064	if (needCleanFake) {
5065	double dummyChangeCost = 0.0;
5066	changeBounds(3, NULL, dummyChangeCost);
5067	}
5068	#if 0
5069	if (objectiveValue_ < lastObjectiveValue_ - 1.0e-8 *
5070	CoinMax(fabs(objectivevalue_), fabs(lastObjectiveValue_))) {
5071	} else {
5072	lastObjectiveValue_ = objectiveValue_;
5073	}
5074	#endif
5075	if (type == 0 \|\| type == 1) {
5076	if (!type && arraysNotCreated) {
5077	// create save arrays
5078	delete [] saveStatus_;
5079	delete [] savedSolution_;
5080	saveStatus_ = new unsigned char [numberRows_+numberColumns_];
5081	savedSolution_ = new double [numberRows_+numberColumns_];
5082	}
5083	// save arrays
5084	CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
5085	CoinMemcpyN(rowActivityWork_,
5086	numberRows_, savedSolution_ + numberColumns_);
5087	CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
5088	// save extra stuff
5089	int dummy;
5090	matrix_->generalExpanded(this, 5, dummy);
5091	}
5092	if (weightsSaved) {
5093	// restore weights (if saved) - also recompute infeasibility list
5094	if (!reallyBadProblems && (largestPrimalError_ < 100.0 \|\| numberPivots > 10)) {
5095	if (tentativeStatus > -3)
5096	dualRowPivot_->saveWeights(this, (type < 2) ? 2 : 4);
5097	else
5098	dualRowPivot_->saveWeights(this, 3);
5099	} else {
5100	// reset weights or scale back
5101	dualRowPivot_->saveWeights(this, 6);
5102	}
5103	} else if (weightsSaved2 && numberPrimalInfeasibilities_) {
5104	dualRowPivot_->saveWeights(this, 3);
5105	}
5106	}
5107	// see if cutoff reached
5108	double limit = 0.0;
5109	getDblParam(ClpDualObjectiveLimit, limit);
5110	#if 0
5111	if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5112	limit + 100.0) {
5113	printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
5114	limit, objectiveValue_, objectiveValue(), computeInternalObjectiveValue(), sumDualInfeasibilities_);
5115	}
5116	#endif
5117	if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5118	limit && !numberAtFakeBound()) {
5119	bool looksInfeasible = !numberDualInfeasibilities_;
5120	if (objectiveValue()optimizationDirection_ > limit + fabs(0.1 limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 &&
5121	sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000)
5122	looksInfeasible = true;
5123	if (looksInfeasible) {
5124	// Even if not perturbed internal costs may have changed
5125	// be careful
5126	if (true \|\| numberIterations_) {
5127	if(computeInternalObjectiveValue() > limit) {
5128	problemStatus_ = 1;
5129	secondaryStatus_ = 1; // and say was on cutoff
5130	}
5131	} else {
5132	problemStatus_ = 1;
5133	secondaryStatus_ = 1; // and say was on cutoff
5134	}
5135	}
5136	}
5137	// If we are in trouble and in branch and bound give up
5138	if ((specialOptions_ & 1024) != 0) {
5139	int looksBad = 0;
5140	if (largestPrimalError_ * largestDualError_ > 1.0e2) {
5141	looksBad = 1;
5142	} else if (largestPrimalError_ > 1.0e-2
5143	&& objectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) {
5144	looksBad = 2;
5145	}
5146	if (looksBad) {
5147	if (factorization_->pivotTolerance() < 0.9) {
5148	// up tolerance
5149	factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
5150	} else if (numberIterations_ > 10000) {
5151	if (handler_->logLevel() > 2)
5152	printf("bad dual - saying infeasible %d\n", looksBad);
5153	problemStatus_ = 1;
5154	secondaryStatus_ = 1; // and say was on cutoff
5155	} else if (largestPrimalError_ > 1.0e5) {
5156	{
5157	int iBigB = -1;
5158	double bigB = 0.0;
5159	int iBigN = -1;
5160	double bigN = 0.0;
5161	for (int i = 0; i < numberRows_ + numberColumns_; i++) {
5162	double value = fabs(solution_[i]);
5163	if (getStatus(i) == basic) {
5164	if (value > bigB) {
5165	bigB = value;
5166	iBigB = i;
5167	}
5168	} else {
5169	if (value > bigN) {
5170	bigN = value;
5171	iBigN = i;
5172	}
5173	}
5174	}
5175	#ifdef CLP_INVESTIGATE
5176	if (bigB > 1.0e8 \|\| bigN > 1.0e8) {
5177	if (handler_->logLevel() > 0)
5178	printf("it %d - basic %d %g, nonbasic %d %g\n",
5179	numberIterations_, iBigB, bigB, iBigN, bigN);
5180	}
5181	#endif
5182	}
5183	#if COIN_DEVELOP!=2
5184	if (handler_->logLevel() > 2)
5185	#endif
5186	printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_);
5187	allSlackBasis(true);
5188	problemStatus_ = 10;
5189	}
5190	}
5191	}
5192	if (problemStatus_ < 0 && !changeMade_) {
5193	problemStatus_ = 4; // unknown
5194	}
5195	lastGoodIteration_ = numberIterations_;
5196	if (numberIterations_ > lastBadIteration_ + 100)
5197	moreSpecialOptions_ &= ~16; // clear check accuracy flag
5198	if (problemStatus_ < 0) {
5199	sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
5200	if (sumDualInfeasibilities_)
5201	numberDualInfeasibilities_ = 1;
5202	}
5203	#ifdef CLP_REPORT_PROGRESS
5204	if (ixxxxxx > ixxyyyy - 3) {
5205	printf("objectiveValue_ %g\n", objectiveValue_);
5206	handler_->setLogLevel(63);
5207	int nTotal = numberColumns_ + numberRows_;
5208	double newObj = 0.0;
5209	for (int i = 0; i < nTotal; i++) {
5210	if (solution_[i])
5211	newObj += solution_[i] * cost_[i];
5212	}
5213	printf("xxx obj %g\n", newObj);
5214	// for now - recompute all
5215	gutsOfSolution(NULL, NULL);
5216	newObj = 0.0;
5217	for (int i = 0; i < nTotal; i++) {
5218	if (solution_[i])
5219	newObj += solution_[i] * cost_[i];
5220	}
5221	printf("yyy obj %g %g\n", newObj, objectiveValue_);
5222	progress_.modifyObjective(objectiveValue_
5223	- bestPossibleImprovement_);
5224	}
5225	#endif
5226	#if 1
5227	double thisObj = progress_.lastObjective(0);
5228	double lastObj = progress_.lastObjective(1);
5229	if (lastObj > thisObj + 1.0e-4 * CoinMax(fabs(thisObj), fabs(lastObj)) + 1.0e-4
5230	&& givenDuals == NULL && firstFree_ < 0) {
5231	int maxFactor = factorization_->maximumPivots();
5232	if (maxFactor > 10) {
5233	if (forceFactorization_ < 0)
5234	forceFactorization_ = maxFactor;
5235	forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
5236	//printf("Reducing factorization frequency\n");
5237	}
5238	}
5239	#endif
5240	// Allow matrices to be sorted etc
5241	int fake = -999; // signal sort
5242	matrix_->correctSequence(this, fake, fake);
5243	if (alphaAccuracy_ > 0.0)
5244	alphaAccuracy_ = 1.0;
5245	// If we are stopping - use plausible objective
5246	// Maybe only in fast dual
5247	if (problemStatus_ > 2)
5248	objectiveValue_ = approximateObjective;
5249	if (problemStatus_ == 1 && (progressFlag_&8) != 0 &&
5250	fabs(objectiveValue_)