Context Navigation

ClpSimplexDual.cpp @ 1785

Last change on this file since 1785 was 1785, checked in by forrest, 9 years ago
patches to try and make parametrics faster
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 326.5 KB

Line
1	/* $Id: ClpSimplexDual.cpp 1785 2011-08-25 10:17:49Z 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	#if CAN_HAVE_ZERO_OBJ>1
1226	if ((specialOptions_&2097152)!=0)
1227	theta_=0.0;
1228	#endif
1229	} else {
1230	// Make sure direction plausible
1231	CoinAssert (upperOut_ < 1.0e50 \|\| lowerOut_ > -1.0e50);
1232	// If in integer cleanup do direction using duals
1233	// may be wrong way round
1234	if(ifValuesPass == 2) {
1235	if (dual_[pivotRow_] > 0.0) {
1236	// this will give a -1 in pivot row (as slacks are -1.0)
1237	directionOut_ = 1;
1238	} else {
1239	directionOut_ = -1;
1240	}
1241	}
1242	if (directionOut_ < 0 && fabs(valueOut_ - upperOut_) > dualBound_ + primalTolerance_) {
1243	if (fabs(valueOut_ - upperOut_) > fabs(valueOut_ - lowerOut_))
1244	directionOut_ = 1;
1245	} else if (directionOut_ > 0 && fabs(valueOut_ - lowerOut_) > dualBound_ + primalTolerance_) {
1246	if (fabs(valueOut_ - upperOut_) < fabs(valueOut_ - lowerOut_))
1247	directionOut_ = -1;
1248	}
1249	double direction = directionOut_;
1250	rowArray_[0]->createPacked(1, &pivotRow_, &direction);
1251	factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
1252	// put row of tableau in rowArray[0] and columnArray[0]
1253	if (!scaledMatrix_) {
1254	matrix_->transposeTimes(this, -1.0,
1255	rowArray_[0], rowArray_[3], columnArray_[0]);
1256	} else {
1257	double * saveR = rowScale_;
1258	double * saveC = columnScale_;
1259	rowScale_ = NULL;
1260	columnScale_ = NULL;
1261	scaledMatrix_->transposeTimes(this, -1.0,
1262	rowArray_[0], rowArray_[3], columnArray_[0]);
1263	rowScale_ = saveR;
1264	columnScale_ = saveC;
1265	}
1266	acceptablePivot *= 10.0;
1267	// do ratio test
1268	if (ifValuesPass == 1) {
1269	checkPossibleValuesMove(rowArray_[0], columnArray_[0],
1270	acceptablePivot);
1271	} else {
1272	checkPossibleCleanup(rowArray_[0], columnArray_[0],
1273	acceptablePivot);
1274	if (sequenceIn_ < 0) {
1275	rowArray_[0]->clear();
1276	columnArray_[0]->clear();
1277	continue; // can't do anything
1278	}
1279	}
1280
1281	// recompute true dualOut_
1282	if (directionOut_ < 0) {
1283	dualOut_ = valueOut_ - upperOut_;
1284	} else {
1285	dualOut_ = lowerOut_ - valueOut_;
1286	}
1287	// check what happened if was values pass
1288	// may want to move part way i.e. movement
1289	bool normalIteration = (sequenceIn_ != sequenceOut_);
1290
1291	clearPivoted(sequenceOut_); // make sure won't be done again
1292	// see if end of values pass
1293	if (!numberCandidates) {
1294	int iRow;
1295	delete [] candidateList;
1296	delete [] givenDuals;
1297	candidate = -2; // -2 signals end
1298	givenDuals = NULL;
1299	candidateList = NULL;
1300	ifValuesPass = 1;
1301	for (iRow = 0; iRow < numberRows_; iRow++) {
1302	int iPivot = pivotVariable_[iRow];
1303	//assert (fabs(dj_[iPivot]),1.0e-5);
1304	clearPivoted(iPivot);
1305	}
1306	}
1307	if (!normalIteration) {
1308	//rowArray_[0]->cleanAndPackSafe(1.0e-60);
1309	//columnArray_[0]->cleanAndPackSafe(1.0e-60);
1310	updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
1311	if (candidate == -2)
1312	problemStatus_ = -2;
1313	continue; // skip rest of iteration
1314	} else {
1315	// recompute dualOut_
1316	if (directionOut_ < 0) {
1317	dualOut_ = valueOut_ - upperOut_;
1318	} else {
1319	dualOut_ = lowerOut_ - valueOut_;
1320	}
1321	}
1322	}
1323	if (sequenceIn_ >= 0) {
1324	// normal iteration
1325	// update the incoming column
1326	double btranAlpha = -alpha_ * directionOut_; // for check
1327	unpackPacked(rowArray_[1]);
1328	// moved into updateWeights - factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
1329	// and update dual weights (can do in parallel - with extra array)
1330	alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
1331	rowArray_[2],
1332	rowArray_[3],
1333	rowArray_[1]);
1334	// see if update stable
1335	#ifdef CLP_DEBUG
1336	if ((handler_->logLevel() & 32))
1337	printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
1338	#endif
1339	double checkValue = 1.0e-7;
1340	// if can't trust much and long way from optimal then relax
1341	if (largestPrimalError_ > 10.0)
1342	checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
1343	if (fabs(btranAlpha) < 1.0e-12 \|\| fabs(alpha_) < 1.0e-12 \|\|
1344	fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
1345	handler_->message(CLP_DUAL_CHECK, messages_)
1346	<< btranAlpha
1347	<< alpha_
1348	<< CoinMessageEol;
1349	if (factorization_->pivots()) {
1350	dualRowPivot_->unrollWeights();
1351	problemStatus_ = -2; // factorize now
1352	rowArray_[0]->clear();
1353	rowArray_[1]->clear();
1354	columnArray_[0]->clear();
1355	returnCode = -2;
1356	break;
1357	} else {
1358	// take on more relaxed criterion
1359	double test;
1360	if (fabs(btranAlpha) < 1.0e-8 \|\| fabs(alpha_) < 1.0e-8)
1361	test = 1.0e-1 * fabs(alpha_);
1362	else
1363	test = 1.0e-4 * (1.0 + fabs(alpha_));
1364	if (fabs(btranAlpha) < 1.0e-12 \|\| fabs(alpha_) < 1.0e-12 \|\|
1365	fabs(btranAlpha - alpha_) > test) {
1366	dualRowPivot_->unrollWeights();
1367	// need to reject something
1368	char x = isColumn(sequenceOut_) ? 'C' : 'R';
1369	handler_->message(CLP_SIMPLEX_FLAG, messages_)
1370	<< x << sequenceWithin(sequenceOut_)
1371	<< CoinMessageEol;
1372	#ifdef COIN_DEVELOP
1373	printf("flag a %g %g\n", btranAlpha, alpha_);
1374	#endif
1375	//#define FEB_TRY
1376	#if 1 //def FEB_TRY
1377	// Make safer?
1378	factorization_->saferTolerances (-0.99, -1.03);
1379	#endif
1380	setFlagged(sequenceOut_);
1381	progress_.clearBadTimes();
1382	lastBadIteration_ = numberIterations_; // say be more cautious
1383	rowArray_[0]->clear();
1384	rowArray_[1]->clear();
1385	columnArray_[0]->clear();
1386	if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
1387	//printf("I think should declare infeasible\n");
1388	problemStatus_ = 1;
1389	returnCode = 1;
1390	break;
1391	}
1392	continue;
1393	}
1394	}
1395	}
1396	// update duals BEFORE replaceColumn so can do updateColumn
1397	double objectiveChange = 0.0;
1398	// do duals first as variables may flip bounds
1399	// rowArray_[0] and columnArray_[0] may have flips
1400	// so use rowArray_[3] for work array from here on
1401	int nswapped = 0;
1402	//rowArray_[0]->cleanAndPackSafe(1.0e-60);
1403	//columnArray_[0]->cleanAndPackSafe(1.0e-60);
1404	if (candidate == -1) {
1405	#if CLP_CAN_HAVE_ZERO_OBJ>1
1406	if ((specialOptions_&2097152)==0) {
1407	#endif
1408	// make sure incoming doesn't count
1409	Status saveStatus = getStatus(sequenceIn_);
1410	setStatus(sequenceIn_, basic);
1411	nswapped = updateDualsInDual(rowArray_[0], columnArray_[0],
1412	rowArray_[2], theta_,
1413	objectiveChange, false);
1414	setStatus(sequenceIn_, saveStatus);
1415	#if CLP_CAN_HAVE_ZERO_OBJ>1
1416	} else {
1417	rowArray_[0]->clear();
1418	rowArray_[2]->clear();
1419	columnArray_[0]->clear();
1420	}
1421	#endif
1422	} else {
1423	updateDualsInValuesPass(rowArray_[0], columnArray_[0], theta_);
1424	}
1425	double oldDualOut = dualOut_;
1426	// which will change basic solution
1427	if (nswapped) {
1428	if (rowArray_[2]->getNumElements()) {
1429	factorization_->updateColumn(rowArray_[3], rowArray_[2]);
1430	dualRowPivot_->updatePrimalSolution(rowArray_[2],
1431	1.0, objectiveChange);
1432	}
1433	// recompute dualOut_
1434	valueOut_ = solution_[sequenceOut_];
1435	if (directionOut_ < 0) {
1436	dualOut_ = valueOut_ - upperOut_;
1437	} else {
1438	dualOut_ = lowerOut_ - valueOut_;
1439	}
1440	#if 0
1441	if (dualOut_ < 0.0) {
1442	#ifdef CLP_DEBUG
1443	if (handler_->logLevel() & 32) {
1444	printf(" dualOut_ %g %g save %g\n", dualOut_, averagePrimalInfeasibility, saveDualOut);
1445	printf("values %g %g %g %g %g %g %g\n", lowerOut_, valueOut_, upperOut_,
1446	objectiveChange,);
1447	}
1448	#endif
1449	if (upperOut_ == lowerOut_)
1450	dualOut_ = 0.0;
1451	}
1452	if(dualOut_ < -CoinMax(1.0e-12 * averagePrimalInfeasibility, 1.0e-8)
1453	&& factorization_->pivots() > 100 &&
1454	getStatus(sequenceIn_) != isFree) {
1455	// going backwards - factorize
1456	dualRowPivot_->unrollWeights();
1457	problemStatus_ = -2; // factorize now
1458	returnCode = -2;
1459	break;
1460	}
1461	#endif
1462	}
1463	// amount primal will move
1464	double movement = -dualOut_ * directionOut_ / alpha_;
1465	double movementOld = oldDualOut * directionOut_ / alpha_;
1466	// so objective should increase by fabs(dj)*movement
1467	// but we already have objective change - so check will be good
1468	if (objectiveChange + fabs(movementOld * dualIn_) < -CoinMax(1.0e-5, 1.0e-12 * fabs(objectiveValue_))) {
1469	#ifdef CLP_DEBUG
1470	if (handler_->logLevel() & 32)
1471	printf("movement %g, swap change %g, rest %g * %g\n",
1472	objectiveChange + fabs(movement * dualIn_),
1473	objectiveChange, movement, dualIn_);
1474	#endif
1475	if(factorization_->pivots()) {
1476	// going backwards - factorize
1477	dualRowPivot_->unrollWeights();
1478	problemStatus_ = -2; // factorize now
1479	returnCode = -2;
1480	break;
1481	}
1482	}
1483	// if stable replace in basis
1484	int updateStatus = factorization_->replaceColumn(this,
1485	rowArray_[2],
1486	rowArray_[1],
1487	pivotRow_,
1488	alpha_,
1489	(moreSpecialOptions_ & 16) != 0,
1490	acceptablePivot);
1491	// If looks like bad pivot - refactorize
1492	if (fabs(dualOut_) > 1.0e50)
1493	updateStatus = 2;
1494	// if no pivots, bad update but reasonable alpha - take and invert
1495	if (updateStatus == 2 &&
1496	!factorization_->pivots() && fabs(alpha_) > 1.0e-5)
1497	updateStatus = 4;
1498	if (updateStatus == 1 \|\| updateStatus == 4) {
1499	// slight error
1500	if (factorization_->pivots() > 5 \|\| updateStatus == 4) {
1501	problemStatus_ = -2; // factorize now
1502	returnCode = -3;
1503	}
1504	} else if (updateStatus == 2) {
1505	// major error
1506	dualRowPivot_->unrollWeights();
1507	// later we may need to unwind more e.g. fake bounds
1508	if (factorization_->pivots() &&
1509	((moreSpecialOptions_ & 16) == 0 \|\| factorization_->pivots() > 4)) {
1510	problemStatus_ = -2; // factorize now
1511	returnCode = -2;
1512	moreSpecialOptions_ \|= 16;
1513	break;
1514	} else {
1515	// need to reject something
1516	char x = isColumn(sequenceOut_) ? 'C' : 'R';
1517	handler_->message(CLP_SIMPLEX_FLAG, messages_)
1518	<< x << sequenceWithin(sequenceOut_)
1519	<< CoinMessageEol;
1520	#ifdef COIN_DEVELOP
1521	printf("flag b %g\n", alpha_);
1522	#endif
1523	setFlagged(sequenceOut_);
1524	progress_.clearBadTimes();
1525	lastBadIteration_ = numberIterations_; // say be more cautious
1526	rowArray_[0]->clear();
1527	rowArray_[1]->clear();
1528	columnArray_[0]->clear();
1529	// make sure dual feasible
1530	// look at all rows and columns
1531	double objectiveChange = 0.0;
1532	updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
1533	0.0, objectiveChange, true);
1534	rowArray_[1]->clear();
1535	columnArray_[0]->clear();
1536	continue;
1537	}
1538	} else if (updateStatus == 3) {
1539	// out of memory
1540	// increase space if not many iterations
1541	if (factorization_->pivots() <
1542	0.5 * factorization_->maximumPivots() &&
1543	factorization_->pivots() < 200)
1544	factorization_->areaFactor(
1545	factorization_->areaFactor() * 1.1);
1546	problemStatus_ = -2; // factorize now
1547	} else if (updateStatus == 5) {
1548	problemStatus_ = -2; // factorize now
1549	}
1550	// update primal solution
1551	if (theta_ < 0.0 && candidate == -1) {
1552	#ifdef CLP_DEBUG
1553	if (handler_->logLevel() & 32)
1554	printf("negative theta %g\n", theta_);
1555	#endif
1556	theta_ = 0.0;
1557	}
1558	// do actual flips
1559	flipBounds(rowArray_[0], columnArray_[0]);
1560	//rowArray_[1]->expand();
1561	dualRowPivot_->updatePrimalSolution(rowArray_[1],
1562	movement,
1563	objectiveChange);
1564	#ifdef CLP_DEBUG
1565	double oldobj = objectiveValue_;
1566	#endif
1567	// modify dualout
1568	dualOut_ /= alpha_;
1569	dualOut_ *= -directionOut_;
1570	//setStatus(sequenceIn_,basic);
1571	dj_[sequenceIn_] = 0.0;
1572	double oldValue = valueIn_;
1573	if (directionIn_ == -1) {
1574	// as if from upper bound
1575	valueIn_ = upperIn_ + dualOut_;
1576	} else {
1577	// as if from lower bound
1578	valueIn_ = lowerIn_ + dualOut_;
1579	}
1580	objectiveChange += cost_[sequenceIn_] * (valueIn_ - oldValue);
1581	// outgoing
1582	// set dj to zero unless values pass
1583	if (directionOut_ > 0) {
1584	valueOut_ = lowerOut_;
1585	if (candidate == -1)
1586	dj_[sequenceOut_] = theta_;
1587	} else {
1588	valueOut_ = upperOut_;
1589	if (candidate == -1)
1590	dj_[sequenceOut_] = -theta_;
1591	}
1592	solution_[sequenceOut_] = valueOut_;
1593	int whatNext = housekeeping(objectiveChange);
1594	#ifdef CLP_REPORT_PROGRESS
1595	if (ixxxxxx > ixxyyyy - 5) {
1596	handler_->setLogLevel(63);
1597	int nTotal = numberColumns_ + numberRows_;
1598	double oldObj = 0.0;
1599	double newObj = 0.0;
1600	for (int i = 0; i < nTotal; i++) {
1601	if (savePSol[i])
1602	oldObj += savePSol[i] * saveCost[i];
1603	if (solution_[i])
1604	newObj += solution_[i] * cost_[i];
1605	bool printIt = false;
1606	if (cost_[i] != saveCost[i])
1607	printIt = true;
1608	if (status_[i] != saveStat[i])
1609	printIt = true;
1610	if (printIt)
1611	printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
1612	i, saveStat[i], saveCost[i], savePSol[i],
1613	status_[i], cost_[i], solution_[i]);
1614	// difference
1615	savePSol[i] = solution_[i] - savePSol[i];
1616	}
1617	printf("pivots %d, old obj %g new %g\n",
1618	factorization_->pivots(),
1619	oldObj, newObj);
1620	memset(saveDj, 0, numberRows_ * sizeof(double));
1621	times(1.0, savePSol, saveDj);
1622	double largest = 1.0e-6;
1623	int k = -1;
1624	for (int i = 0; i < numberRows_; i++) {
1625	saveDj[i] -= savePSol[i+numberColumns_];
1626	if (fabs(saveDj[i]) > largest) {
1627	largest = fabs(saveDj[i]);
1628	k = i;
1629	}
1630	}
1631	if (k >= 0)
1632	printf("Not null %d %g\n", k, largest);
1633	}
1634	#endif
1635	#ifdef VUB
1636	{
1637	if ((sequenceIn_ < numberColumns_ && vub[sequenceIn_] >= 0) \|\| toVub[sequenceIn_] >= 0 \|\|
1638	(sequenceOut_ < numberColumns_ && vub[sequenceOut_] >= 0) \|\| toVub[sequenceOut_] >= 0) {
1639	int inSequence = sequenceIn_;
1640	int inVub = -1;
1641	if (sequenceIn_ < numberColumns_)
1642	inVub = vub[sequenceIn_];
1643	int inBack = toVub[inSequence];
1644	int inSlack = -1;
1645	if (inSequence >= numberColumns_ && inBack >= 0) {
1646	inSlack = inSequence - numberColumns_;
1647	inSequence = inBack;
1648	inBack = toVub[inSequence];
1649	}
1650	if (inVub >= 0)
1651	printf("Vub %d in ", inSequence);
1652	if (inBack >= 0 && inSlack < 0)
1653	printf("%d (descendent of %d) in ", inSequence, inBack);
1654	if (inSlack >= 0)
1655	printf("slack for row %d -> %d (descendent of %d) in ", inSlack, inSequence, inBack);
1656	int outSequence = sequenceOut_;
1657	int outVub = -1;
1658	if (sequenceOut_ < numberColumns_)
1659	outVub = vub[sequenceOut_];
1660	int outBack = toVub[outSequence];
1661	int outSlack = -1;
1662	if (outSequence >= numberColumns_ && outBack >= 0) {
1663	outSlack = outSequence - numberColumns_;
1664	outSequence = outBack;
1665	outBack = toVub[outSequence];
1666	}
1667	if (outVub >= 0)
1668	printf("Vub %d out ", outSequence);
1669	if (outBack >= 0 && outSlack < 0)
1670	printf("%d (descendent of %d) out ", outSequence, outBack);
1671	if (outSlack >= 0)
1672	printf("slack for row %d -> %d (descendent of %d) out ", outSlack, outSequence, outBack);
1673	printf("\n");
1674	}
1675	}
1676	#endif
1677	#if 0
1678	if (numberIterations_ > 206033)
1679	handler_->setLogLevel(63);
1680	if (numberIterations_ > 210567)
1681	exit(77);
1682	#endif
1683	if (!givenDuals && ifValuesPass && ifValuesPass != 2) {
1684	handler_->message(CLP_END_VALUES_PASS, messages_)
1685	<< numberIterations_;
1686	whatNext = 1;
1687	}
1688	#ifdef CHECK_ACCURACY
1689	if (whatNext) {
1690	CoinMemcpyN(solution_, (numberRows_ + numberColumns_), zzzzzz);
1691	}
1692	#endif
1693	//if (numberIterations_==1890)
1694	//whatNext=1;
1695	//if (numberIterations_>2000)
1696	//exit(77);
1697	// and set bounds correctly
1698	originalBound(sequenceIn_);
1699	changeBound(sequenceOut_);
1700	#ifdef CLP_DEBUG
1701	if (objectiveValue_ < oldobj - 1.0e-5 && (handler_->logLevel() & 16))
1702	printf("obj backwards %g %g\n", objectiveValue_, oldobj);
1703	#endif
1704	#if 0
1705	{
1706	for (int i = 0; i < numberRows_ + numberColumns_; i++) {
1707	FakeBound bound = getFakeBound(i);
1708	if (bound == ClpSimplexDual::upperFake) {
1709	assert (upper_[i] < 1.0e20);
1710	} else if (bound == ClpSimplexDual::lowerFake) {
1711	assert (lower_[i] > -1.0e20);
1712	} else if (bound == ClpSimplexDual::bothFake) {
1713	assert (upper_[i] < 1.0e20);
1714	assert (lower_[i] > -1.0e20);
1715	}
1716	}
1717	}
1718	#endif
1719	if (whatNext == 1 \|\| candidate == -2) {
1720	problemStatus_ = -2; // refactorize
1721	} else if (whatNext == 2) {
1722	// maximum iterations or equivalent
1723	problemStatus_ = 3;
1724	returnCode = 3;
1725	break;
1726	}
1727	// Check event
1728	{
1729	int status = eventHandler_->event(ClpEventHandler::endOfIteration);
1730	if (status >= 0) {
1731	problemStatus_ = 5;
1732	secondaryStatus_ = ClpEventHandler::endOfIteration;
1733	returnCode = 4;
1734	break;
1735	}
1736	}
1737	} else {
1738	#ifdef CLP_INVESTIGATE_SERIAL
1739	z_thinks = 1;
1740	#endif
1741	// no incoming column is valid
1742	pivotRow_ = -1;
1743	#ifdef CLP_DEBUG
1744	if (handler_->logLevel() & 32)
1745	printf("** no column pivot\n");
1746	#endif
1747	if (factorization_->pivots() < 2 && acceptablePivot_ <= 1.0e-8) {
1748	//&&goodAccuracy()) {
1749	// If not in branch and bound etc save ray
1750	delete [] ray_;
1751	if ((specialOptions_&(1024 \| 4096)) == 0 \|\| (specialOptions_ & 32) != 0) {
1752	// create ray anyway
1753	ray_ = new double [ numberRows_];
1754	rowArray_[0]->expand(); // in case packed
1755	CoinMemcpyN(rowArray_[0]->denseVector(), numberRows_, ray_);
1756	} else {
1757	ray_ = NULL;
1758	}
1759	// If we have just factorized and infeasibility reasonable say infeas
1760	double dualTest = ((specialOptions_ & 4096) != 0) ? 1.0e8 : 1.0e13;
1761	if (((specialOptions_ & 4096) != 0 \|\| bestPossiblePivot < 1.0e-11) && dualBound_ > dualTest) {
1762	double testValue = 1.0e-4;
1763	if (!factorization_->pivots() && numberPrimalInfeasibilities_ == 1)
1764	testValue = 1.0e-6;
1765	if (valueOut_ > upperOut_ + testValue \|\| valueOut_ < lowerOut_ - testValue
1766	\|\| (specialOptions_ & 64) == 0) {
1767	// say infeasible
1768	problemStatus_ = 1;
1769	// unless primal feasible!!!!
1770	//printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
1771	// numberDualInfeasibilities_,sumDualInfeasibilities_);
1772	//#define TEST_CLP_NODE
1773	#ifndef TEST_CLP_NODE
1774	// Should be correct - but ...
1775	int numberFake = numberAtFakeBound();
1776	double sumPrimal = (!numberFake) ? 2.0e5 : sumPrimalInfeasibilities_;
1777	if (sumPrimalInfeasibilities_ < 1.0e-3 \|\| sumDualInfeasibilities_ > 1.0e-5 \|\|
1778	(sumPrimal < 1.0e5 && (specialOptions_ & 1024) != 0 && factorization_->pivots())) {
1779	if (sumPrimal > 50.0 && factorization_->pivots() > 2) {
1780	problemStatus_ = -4;
1781	#ifdef COIN_DEVELOP
1782	printf("status to -4 at %d - primalinf %g pivots %d\n",
1783	__LINE__, sumPrimalInfeasibilities_,
1784	factorization_->pivots());
1785	#endif
1786	} else {
1787	problemStatus_ = 10;
1788	#if COIN_DEVELOP>1
1789	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",
1790	__LINE__, numberPrimalInfeasibilities_,
1791	sumPrimalInfeasibilities_,
1792	numberDualInfeasibilities_, sumDualInfeasibilities_,
1793	numberFake_, dualBound_, factorization_->pivots(),
1794	(specialOptions_ & 1024) != 0 ? 1 : 0,
1795	(specialOptions_ & 2048) != 0 ? 1 : 0,
1796	(specialOptions_ & 4096) != 0 ? 1 : 0,
1797	(specialOptions_ & 8192) != 0 ? 1 : 0,
1798	(specialOptions_ & 16384) != 0 ? 1 : 0
1799	);
1800	#endif
1801	// Get rid of objective
1802	if ((specialOptions_ & 16384) == 0)
1803	objective_ = new ClpLinearObjective(NULL, numberColumns_);
1804	}
1805	}
1806	#else
1807	if (sumPrimalInfeasibilities_ < 1.0e-3 \|\| sumDualInfeasibilities_ > 1.0e-6) {
1808	#ifdef COIN_DEVELOP
1809	printf("at %d - primal %d %g - dual %d %g fake %d weight %g - pivs %d\n",
1810	__LINE__, numberPrimalInfeasibilities_,
1811	sumPrimalInfeasibilities_,
1812	numberDualInfeasibilities_, sumDualInfeasibilities_,
1813	numberFake_, dualBound_, factorization_->pivots());
1814	#endif
1815	if ((specialOptions_ & 1024) != 0 && factorization_->pivots()) {
1816	problemStatus_ = 10;
1817	#if COIN_DEVELOP>1
1818	printf("returning at %d\n", __LINE__);
1819	#endif
1820	// Get rid of objective
1821	if ((specialOptions_ & 16384) == 0)
1822	objective_ = new ClpLinearObjective(NULL, numberColumns_);
1823	}
1824	}
1825	#endif
1826	rowArray_[0]->clear();
1827	columnArray_[0]->clear();
1828	returnCode = 1;
1829	break;
1830	}
1831	}
1832	// If special option set - put off as long as possible
1833	if ((specialOptions_ & 64) == 0 \|\| (moreSpecialOptions_ & 64) != 0) {
1834	if (factorization_->pivots() == 0)
1835	problemStatus_ = -4; //say looks infeasible
1836	} else {
1837	// flag
1838	char x = isColumn(sequenceOut_) ? 'C' : 'R';
1839	handler_->message(CLP_SIMPLEX_FLAG, messages_)
1840	<< x << sequenceWithin(sequenceOut_)
1841	<< CoinMessageEol;
1842	#ifdef COIN_DEVELOP
1843	printf("flag c\n");
1844	#endif
1845	setFlagged(sequenceOut_);
1846	if (!factorization_->pivots()) {
1847	rowArray_[0]->clear();
1848	columnArray_[0]->clear();
1849	continue;
1850	}
1851	}
1852	}
1853	if (factorization_->pivots() < 5 && acceptablePivot_ > 1.0e-8)
1854	acceptablePivot_ = 1.0e-8;
1855	rowArray_[0]->clear();
1856	columnArray_[0]->clear();
1857	returnCode = 1;
1858	break;
1859	}
1860	} else {
1861	#ifdef CLP_INVESTIGATE_SERIAL
1862	z_thinks = 0;
1863	#endif
1864	// no pivot row
1865	#ifdef CLP_DEBUG
1866	if (handler_->logLevel() & 32)
1867	printf("** no row pivot\n");
1868	#endif
1869	// If in branch and bound try and get rid of fixed variables
1870	if ((specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
1871	assert (!candidateList);
1872	candidateList = new int[numberRows_];
1873	int iRow;
1874	for (iRow = 0; iRow < numberRows_; iRow++) {
1875	int iPivot = pivotVariable_[iRow];
1876	if (flagged(iPivot) \|\| !pivoted(iPivot))
1877	continue;
1878	assert (iPivot < numberColumns_ && lower_[iPivot] == upper_[iPivot]);
1879	candidateList[numberCandidates++] = iRow;
1880	}
1881	// and set first candidate
1882	if (!numberCandidates) {
1883	delete [] candidateList;
1884	candidateList = NULL;
1885	int iRow;
1886	for (iRow = 0; iRow < numberRows_; iRow++) {
1887	int iPivot = pivotVariable_[iRow];
1888	clearPivoted(iPivot);
1889	}
1890	} else {
1891	ifValuesPass = 2;
1892	continue;
1893	}
1894	}
1895	int numberPivots = factorization_->pivots();
1896	bool specialCase;
1897	int useNumberFake;
1898	returnCode = 0;
1899	if (numberPivots <= CoinMax(dontFactorizePivots_, 20) &&
1900	(specialOptions_ & 2048) != 0 && (true \|\| !numberChanged_ \|\| perturbation_ == 101)
1901	&& dualBound_ >= 1.0e8) {
1902	specialCase = true;
1903	// as dual bound high - should be okay
1904	useNumberFake = 0;
1905	} else {
1906	specialCase = false;
1907	useNumberFake = numberFake_;
1908	}
1909	if (!numberPivots \|\| specialCase) {
1910	// may have crept through - so may be optimal
1911	// check any flagged variables
1912	int iRow;
1913	for (iRow = 0; iRow < numberRows_; iRow++) {
1914	int iPivot = pivotVariable_[iRow];
1915	if (flagged(iPivot))
1916	break;
1917	}
1918	if (iRow < numberRows_ && numberPivots) {
1919	// try factorization
1920	returnCode = -2;
1921	}
1922
1923	if (useNumberFake \|\| numberDualInfeasibilities_) {
1924	// may be dual infeasible
1925	if ((specialOptions_ & 1024) == 0)
1926	problemStatus_ = -5;
1927	else if (!useNumberFake && numberPrimalInfeasibilities_
1928	&& !numberPivots)
1929	problemStatus_ = 1;
1930	} else {
1931	if (iRow < numberRows_) {
1932	#ifdef COIN_DEVELOP
1933	std::cout << "Flagged variables at end - infeasible?" << std::endl;
1934	printf("Probably infeasible - pivot was %g\n", alpha_);
1935	#endif
1936	//if (fabs(alpha_)<1.0e-4) {
1937	//problemStatus_=1;
1938	//} else {
1939	#ifdef CLP_DEBUG
1940	abort();
1941	#endif
1942	//}
1943	problemStatus_ = -5;
1944	} else {
1945	problemStatus_ = 0;
1946	#ifndef CLP_CHECK_NUMBER_PIVOTS
1947	#define CLP_CHECK_NUMBER_PIVOTS 10
1948	#endif
1949	#if CLP_CHECK_NUMBER_PIVOTS < 20
1950	if (numberPivots > CLP_CHECK_NUMBER_PIVOTS) {
1951	#ifndef NDEBUG_CLP
1952	int nTotal = numberRows_ + numberColumns_;
1953	double * comp = CoinCopyOfArray(solution_, nTotal);
1954	#endif
1955	computePrimals(rowActivityWork_, columnActivityWork_);
1956	#ifndef NDEBUG_CLP
1957	double largest = 1.0e-5;
1958	int bad = -1;
1959	for (int i = 0; i < nTotal; i++) {
1960	double value = solution_[i];
1961	double larger = CoinMax(fabs(value), fabs(comp[i]));
1962	double tol = 1.0e-5 + 1.0e-5 * larger;
1963	double diff = fabs(value - comp[i]);
1964	if (diff - tol > largest) {
1965	bad = i;
1966	largest = diff - tol;
1967	}
1968	}
1969	if (bad >= 0)
1970	COIN_DETAIL_PRINT(printf("bad %d old %g new %g\n", bad, comp[bad], solution_[bad]));
1971	#endif
1972	checkPrimalSolution(rowActivityWork_, columnActivityWork_);
1973	if (numberPrimalInfeasibilities_) {
1974	#ifdef CLP_INVESTIGATE
1975	printf("XXX Infeas ? %d inf summing to %g\n", numberPrimalInfeasibilities_,
1976	sumPrimalInfeasibilities_);
1977	#endif
1978	problemStatus_ = -1;
1979	returnCode = -2;
1980	}
1981	#ifndef NDEBUG_CLP
1982	memcpy(solution_, comp, nTotal * sizeof(double));
1983	delete [] comp;
1984	#endif
1985	}
1986	#endif
1987	if (!problemStatus_) {
1988	// make it look OK
1989	numberPrimalInfeasibilities_ = 0;
1990	sumPrimalInfeasibilities_ = 0.0;
1991	numberDualInfeasibilities_ = 0;
1992	sumDualInfeasibilities_ = 0.0;
1993	// May be perturbed
1994	if (perturbation_ == 101 \|\| numberChanged_) {
1995	numberChanged_ = 0; // Number of variables with changed costs
1996	perturbation_ = 102; // stop any perturbations
1997	//double changeCost;
1998	//changeBounds(1,NULL,changeCost);
1999	createRim4(false);
2000	// make sure duals are current
2001	computeDuals(givenDuals);
2002	checkDualSolution();
2003	progress_.modifyObjective(-COIN_DBL_MAX);
2004	if (numberDualInfeasibilities_) {
2005	problemStatus_ = 10; // was -3;
2006	} else {
2007	computeObjectiveValue(true);
2008	}
2009	} else if (numberPivots) {
2010	computeObjectiveValue(true);
2011	}
2012	if (numberPivots < -1000) {
2013	// objective may be wrong
2014	objectiveValue_ = innerProduct(cost_, numberColumns_ + numberRows_, solution_);
2015	objectiveValue_ += objective_->nonlinearOffset();
2016	objectiveValue_ /= (objectiveScale_ * rhsScale_);
2017	if ((specialOptions_ & 16384) == 0) {
2018	// and dual_ may be wrong (i.e. for fixed or basic)
2019	CoinIndexedVector * arrayVector = rowArray_[1];
2020	arrayVector->clear();
2021	int iRow;
2022	double * array = arrayVector->denseVector();
2023	/* Use dual_ instead of array
2024	Even though dual_ is only numberRows_ long this is
2025	okay as gets permuted to longer rowArray_[2]
2026	*/
2027	arrayVector->setDenseVector(dual_);
2028	int * index = arrayVector->getIndices();
2029	int number = 0;
2030	for (iRow = 0; iRow < numberRows_; iRow++) {
2031	int iPivot = pivotVariable_[iRow];
2032	double value = cost_[iPivot];
2033	dual_[iRow] = value;
2034	if (value) {
2035	index[number++] = iRow;
2036	}
2037	}
2038	arrayVector->setNumElements(number);
2039	// Extended duals before "updateTranspose"
2040	matrix_->dualExpanded(this, arrayVector, NULL, 0);
2041	// Btran basic costs
2042	rowArray_[2]->clear();
2043	factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
2044	// and return vector
2045	arrayVector->setDenseVector(array);
2046	}
2047	}
2048	sumPrimalInfeasibilities_ = 0.0;
2049	}
2050	if ((specialOptions_&(1024 + 16384)) != 0 && !problemStatus_) {
2051	CoinIndexedVector * arrayVector = rowArray_[1];
2052	arrayVector->clear();
2053	double * rhs = arrayVector->denseVector();
2054	times(1.0, solution_, rhs);
2055	#ifdef CHECK_ACCURACY
2056	bool bad = false;
2057	#endif
2058	bool bad2 = false;
2059	int i;
2060	for ( i = 0; i < numberRows_; i++) {
2061	if (rhs[i] < rowLowerWork_[i] - primalTolerance_ \|\|
2062	rhs[i] > rowUpperWork_[i] + primalTolerance_) {
2063	bad2 = true;
2064	#ifdef CHECK_ACCURACY
2065	printf("row %d out of bounds %g, %g correct %g bad %g\n", i,
2066	rowLowerWork_[i], rowUpperWork_[i],
2067	rhs[i], rowActivityWork_[i]);
2068	#endif
2069	} else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
2070	#ifdef CHECK_ACCURACY
2071	bad = true;
2072	printf("row %d correct %g bad %g\n", i, rhs[i], rowActivityWork_[i]);
2073	#endif
2074	}
2075	rhs[i] = 0.0;
2076	}
2077	for ( i = 0; i < numberColumns_; i++) {
2078	if (solution_[i] < columnLowerWork_[i] - primalTolerance_ \|\|
2079	solution_[i] > columnUpperWork_[i] + primalTolerance_) {
2080	bad2 = true;
2081	#ifdef CHECK_ACCURACY
2082	printf("column %d out of bounds %g, %g correct %g bad %g\n", i,
2083	columnLowerWork_[i], columnUpperWork_[i],
2084	solution_[i], columnActivityWork_[i]);
2085	#endif
2086	}
2087	}
2088	if (bad2) {
2089	problemStatus_ = -3;
2090	returnCode = -2;
2091	// Force to re-factorize early next time
2092	int numberPivots = factorization_->pivots();
2093	forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2094	}
2095	}
2096	}
2097	}
2098	} else {
2099	problemStatus_ = -3;
2100	returnCode = -2;
2101	// Force to re-factorize early next time
2102	int numberPivots = factorization_->pivots();
2103	forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
2104	}
2105	break;
2106	}
2107	}
2108	if (givenDuals) {
2109	CoinMemcpyN(dj_, numberRows_ + numberColumns_, givenDuals);
2110	// get rid of any values pass array
2111	delete [] candidateList;
2112	}
2113	delete [] dubiousWeights;
2114	#ifdef CLP_REPORT_PROGRESS
2115	if (ixxxxxx > ixxyyyy - 5) {
2116	int nTotal = numberColumns_ + numberRows_;
2117	double oldObj = 0.0;
2118	double newObj = 0.0;
2119	for (int i = 0; i < nTotal; i++) {
2120	if (savePSol[i])
2121	oldObj += savePSol[i] * saveCost[i];
2122	if (solution_[i])
2123	newObj += solution_[i] * cost_[i];
2124	bool printIt = false;
2125	if (cost_[i] != saveCost[i])
2126	printIt = true;
2127	if (status_[i] != saveStat[i])
2128	printIt = true;
2129	if (printIt)
2130	printf("%d old %d cost %g sol %g, new %d cost %g sol %g\n",
2131	i, saveStat[i], saveCost[i], savePSol[i],
2132	status_[i], cost_[i], solution_[i]);
2133	// difference
2134	savePSol[i] = solution_[i] - savePSol[i];
2135	}
2136	printf("exit pivots %d, old obj %g new %g\n",
2137	factorization_->pivots(),
2138	oldObj, newObj);
2139	memset(saveDj, 0, numberRows_ * sizeof(double));
2140	times(1.0, savePSol, saveDj);
2141	double largest = 1.0e-6;
2142	int k = -1;
2143	for (int i = 0; i < numberRows_; i++) {
2144	saveDj[i] -= savePSol[i+numberColumns_];
2145	if (fabs(saveDj[i]) > largest) {
2146	largest = fabs(saveDj[i]);
2147	k = i;
2148	}
2149	}
2150	if (k >= 0)
2151	printf("Not null %d %g\n", k, largest);
2152	}
2153	delete [] savePSol ;
2154	delete [] saveDj ;
2155	delete [] saveCost ;
2156	delete [] saveStat ;
2157	#endif
2158	return returnCode;
2159	}
2160	/* The duals are updated by the given arrays.
2161	Returns number of infeasibilities.
2162	rowArray and columnarray will have flipped
2163	The output vector has movement (row length array) */
2164	int
2165	ClpSimplexDual::updateDualsInDual(CoinIndexedVector * rowArray,
2166	CoinIndexedVector * columnArray,
2167	CoinIndexedVector * outputArray,
2168	double theta,
2169	double & objectiveChange,
2170	bool fullRecompute)
2171	{
2172
2173	outputArray->clear();
2174
2175
2176	int numberInfeasibilities = 0;
2177	int numberRowInfeasibilities = 0;
2178
2179	// get a tolerance
2180	double tolerance = dualTolerance_;
2181	// we can't really trust infeasibilities if there is dual error
2182	double error = CoinMin(1.0e-2, largestDualError_);
2183	// allow tolerance at least slightly bigger than standard
2184	tolerance = tolerance + error;
2185
2186	double changeObj = 0.0;
2187
2188	// Coding is very similar but we can save a bit by splitting
2189	// Do rows
2190	if (!fullRecompute) {
2191	int i;
2192	double * COIN_RESTRICT reducedCost = djRegion(0);
2193	const double * COIN_RESTRICT lower = lowerRegion(0);
2194	const double * COIN_RESTRICT upper = upperRegion(0);
2195	const double * COIN_RESTRICT cost = costRegion(0);
2196	double * COIN_RESTRICT work;
2197	int number;
2198	int * COIN_RESTRICT which;
2199	const unsigned char * COIN_RESTRICT statusArray = status_ + numberColumns_;
2200	assert(rowArray->packedMode());
2201	work = rowArray->denseVector();
2202	number = rowArray->getNumElements();
2203	which = rowArray->getIndices();
2204	double multiplier[] = { -1.0, 1.0};
2205	for (i = 0; i < number; i++) {
2206	int iSequence = which[i];
2207	double alphaI = work[i];
2208	work[i] = 0.0;
2209	int iStatus = (statusArray[iSequence] & 3) - 1;
2210	if (iStatus) {
2211	double value = reducedCost[iSequence] - theta * alphaI;
2212	reducedCost[iSequence] = value;
2213	double mult = multiplier[iStatus-1];
2214	value *= mult;
2215	if (value < -tolerance) {
2216	// flipping bounds
2217	double movement = mult * (lower[iSequence] - upper[iSequence]);
2218	which[numberInfeasibilities++] = iSequence;
2219	#ifndef NDEBUG
2220	if (fabs(movement) >= 1.0e30)
2221	resetFakeBounds(-1000 - iSequence);
2222	#endif
2223	#ifdef CLP_DEBUG
2224	if ((handler_->logLevel() & 32))
2225	printf("%d %d, new dj %g, alpha %g, movement %g\n",
2226	0, iSequence, value, alphaI, movement);
2227	#endif
2228	changeObj -= movement * cost[iSequence];
2229	outputArray->quickAdd(iSequence, movement);
2230	}
2231	}
2232	}
2233	// Do columns
2234	reducedCost = djRegion(1);
2235	lower = lowerRegion(1);
2236	upper = upperRegion(1);
2237	cost = costRegion(1);
2238	// set number of infeasibilities in row array
2239	numberRowInfeasibilities = numberInfeasibilities;
2240	rowArray->setNumElements(numberInfeasibilities);
2241	numberInfeasibilities = 0;
2242	work = columnArray->denseVector();
2243	number = columnArray->getNumElements();
2244	which = columnArray->getIndices();
2245	if ((moreSpecialOptions_ & 8) != 0) {
2246	const unsigned char * COIN_RESTRICT statusArray = status_;
2247	for (i = 0; i < number; i++) {
2248	int iSequence = which[i];
2249	double alphaI = work[i];
2250	work[i] = 0.0;
2251
2252	int iStatus = (statusArray[iSequence] & 3) - 1;
2253	if (iStatus) {
2254	double value = reducedCost[iSequence] - theta * alphaI;
2255	reducedCost[iSequence] = value;
2256	double mult = multiplier[iStatus-1];
2257	value *= mult;
2258	if (value < -tolerance) {
2259	// flipping bounds
2260	double movement = mult * (upper[iSequence] - lower[iSequence]);
2261	which[numberInfeasibilities++] = iSequence;
2262	#ifndef NDEBUG
2263	if (fabs(movement) >= 1.0e30)
2264	resetFakeBounds(-1000 - iSequence);
2265	#endif
2266	#ifdef CLP_DEBUG
2267	if ((handler_->logLevel() & 32))
2268	printf("%d %d, new dj %g, alpha %g, movement %g\n",
2269	1, iSequence, value, alphaI, movement);
2270	#endif
2271	changeObj += movement * cost[iSequence];
2272	matrix_->add(this, outputArray, iSequence, movement);
2273	}
2274	}
2275	}
2276	} else {
2277	for (i = 0; i < number; i++) {
2278	int iSequence = which[i];
2279	double alphaI = work[i];
2280	work[i] = 0.0;
2281
2282	Status status = getStatus(iSequence);
2283	if (status == atLowerBound) {
2284	double value = reducedCost[iSequence] - theta * alphaI;
2285	reducedCost[iSequence] = value;
2286	double movement = 0.0;
2287
2288	if (value < -tolerance) {
2289	// to upper bound
2290	which[numberInfeasibilities++] = iSequence;
2291	movement = upper[iSequence] - lower[iSequence];
2292	#ifndef NDEBUG
2293	if (fabs(movement) >= 1.0e30)
2294	resetFakeBounds(-1000 - iSequence);
2295	#endif
2296	#ifdef CLP_DEBUG
2297	if ((handler_->logLevel() & 32))
2298	printf("%d %d, new dj %g, alpha %g, movement %g\n",
2299	1, iSequence, value, alphaI, movement);
2300	#endif
2301	changeObj += movement * cost[iSequence];
2302	matrix_->add(this, outputArray, iSequence, movement);
2303	}
2304	} else if (status == atUpperBound) {
2305	double value = reducedCost[iSequence] - theta * alphaI;
2306	reducedCost[iSequence] = value;
2307	double movement = 0.0;
2308
2309	if (value > tolerance) {
2310	// to lower bound (if swap)
2311	which[numberInfeasibilities++] = iSequence;
2312	movement = lower[iSequence] - upper[iSequence];
2313	#ifndef NDEBUG
2314	if (fabs(movement) >= 1.0e30)
2315	resetFakeBounds(-1000 - iSequence);
2316	#endif
2317	#ifdef CLP_DEBUG
2318	if ((handler_->logLevel() & 32))
2319	printf("%d %d, new dj %g, alpha %g, movement %g\n",
2320	1, iSequence, value, alphaI, movement);
2321	#endif
2322	changeObj += movement * cost[iSequence];
2323	matrix_->add(this, outputArray, iSequence, movement);
2324	}
2325	} else if (status == isFree) {
2326	double value = reducedCost[iSequence] - theta * alphaI;
2327	reducedCost[iSequence] = value;
2328	}
2329	}
2330	}
2331	} else {
2332	double * COIN_RESTRICT solution = solutionRegion(0);
2333	double * COIN_RESTRICT reducedCost = djRegion(0);
2334	const double * COIN_RESTRICT lower = lowerRegion(0);
2335	const double * COIN_RESTRICT upper = upperRegion(0);
2336	const double * COIN_RESTRICT cost = costRegion(0);
2337	int * COIN_RESTRICT which;
2338	which = rowArray->getIndices();
2339	int iSequence;
2340	for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2341	double value = reducedCost[iSequence];
2342
2343	Status status = getStatus(iSequence + numberColumns_);
2344	// more likely to be at upper bound ?
2345	if (status == atUpperBound) {
2346	double movement = 0.0;
2347	//#define NO_SWAP7
2348	if (value > tolerance) {
2349	// to lower bound (if swap)
2350	// put back alpha
2351	which[numberInfeasibilities++] = iSequence;
2352	movement = lower[iSequence] - upper[iSequence];
2353	changeObj += movement * cost[iSequence];
2354	outputArray->quickAdd(iSequence, -movement);
2355	#ifndef NDEBUG
2356	if (fabs(movement) >= 1.0e30)
2357	resetFakeBounds(-1000 - iSequence);
2358	#endif
2359	#ifndef NO_SWAP7
2360	} else if (value > -tolerance) {
2361	// at correct bound but may swap
2362	FakeBound bound = getFakeBound(iSequence + numberColumns_);
2363	if (bound == ClpSimplexDual::upperFake) {
2364	movement = lower[iSequence] - upper[iSequence];
2365	#ifndef NDEBUG
2366	if (fabs(movement) >= 1.0e30)
2367	resetFakeBounds(-1000 - iSequence);
2368	#endif
2369	setStatus(iSequence + numberColumns_, atLowerBound);
2370	solution[iSequence] = lower[iSequence];
2371	changeObj += movement * cost[iSequence];
2372	//numberFake_--;
2373	//setFakeBound(iSequence+numberColumns_,noFake);
2374	}
2375	#endif
2376	}
2377	} else if (status == atLowerBound) {
2378	double movement = 0.0;
2379
2380	if (value < -tolerance) {
2381	// to upper bound
2382	// put back alpha
2383	which[numberInfeasibilities++] = iSequence;
2384	movement = upper[iSequence] - lower[iSequence];
2385	#ifndef NDEBUG
2386	if (fabs(movement) >= 1.0e30)
2387	resetFakeBounds(-1000 - iSequence);
2388	#endif
2389	changeObj += movement * cost[iSequence];
2390	outputArray->quickAdd(iSequence, -movement);
2391	#ifndef NO_SWAP7
2392	} else if (value < tolerance) {
2393	// at correct bound but may swap
2394	FakeBound bound = getFakeBound(iSequence + numberColumns_);
2395	if (bound == ClpSimplexDual::lowerFake) {
2396	movement = upper[iSequence] - lower[iSequence];
2397	#ifndef NDEBUG
2398	if (fabs(movement) >= 1.0e30)
2399	resetFakeBounds(-1000 - iSequence);
2400	#endif
2401	setStatus(iSequence + numberColumns_, atUpperBound);
2402	solution[iSequence] = upper[iSequence];
2403	changeObj += movement * cost[iSequence];
2404	//numberFake_--;
2405	//setFakeBound(iSequence+numberColumns_,noFake);
2406	}
2407	#endif
2408	}
2409	}
2410	}
2411	// Do columns
2412	solution = solutionRegion(1);
2413	reducedCost = djRegion(1);
2414	lower = lowerRegion(1);
2415	upper = upperRegion(1);
2416	cost = costRegion(1);
2417	// set number of infeasibilities in row array
2418	numberRowInfeasibilities = numberInfeasibilities;
2419	rowArray->setNumElements(numberInfeasibilities);
2420	numberInfeasibilities = 0;
2421	which = columnArray->getIndices();
2422	for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2423	double value = reducedCost[iSequence];
2424
2425	Status status = getStatus(iSequence);
2426	if (status == atLowerBound) {
2427	double movement = 0.0;
2428
2429	if (value < -tolerance) {
2430	// to upper bound
2431	// put back alpha
2432	which[numberInfeasibilities++] = iSequence;
2433	movement = upper[iSequence] - lower[iSequence];
2434	#ifndef NDEBUG
2435	if (fabs(movement) >= 1.0e30)
2436	resetFakeBounds(-1000 - iSequence);
2437	#endif
2438	changeObj += movement * cost[iSequence];
2439	matrix_->add(this, outputArray, iSequence, movement);
2440	#ifndef NO_SWAP7
2441	} else if (value < tolerance) {
2442	// at correct bound but may swap
2443	FakeBound bound = getFakeBound(iSequence);
2444	if (bound == ClpSimplexDual::lowerFake) {
2445	movement = upper[iSequence] - lower[iSequence];
2446	#ifndef NDEBUG
2447	if (fabs(movement) >= 1.0e30)
2448	resetFakeBounds(-1000 - iSequence);
2449	#endif
2450	setStatus(iSequence, atUpperBound);
2451	solution[iSequence] = upper[iSequence];
2452	changeObj += movement * cost[iSequence];
2453	//numberFake_--;
2454	//setFakeBound(iSequence,noFake);
2455	}
2456	#endif
2457	}
2458	} else if (status == atUpperBound) {
2459	double movement = 0.0;
2460
2461	if (value > tolerance) {
2462	// to lower bound (if swap)
2463	// put back alpha
2464	which[numberInfeasibilities++] = iSequence;
2465	movement = lower[iSequence] - upper[iSequence];
2466	#ifndef NDEBUG
2467	if (fabs(movement) >= 1.0e30)
2468	resetFakeBounds(-1000 - iSequence);
2469	#endif
2470	changeObj += movement * cost[iSequence];
2471	matrix_->add(this, outputArray, iSequence, movement);
2472	#ifndef NO_SWAP7
2473	} else if (value > -tolerance) {
2474	// at correct bound but may swap
2475	FakeBound bound = getFakeBound(iSequence);
2476	if (bound == ClpSimplexDual::upperFake) {
2477	movement = lower[iSequence] - upper[iSequence];
2478	#ifndef NDEBUG
2479	if (fabs(movement) >= 1.0e30)
2480	resetFakeBounds(-1000 - iSequence);
2481	#endif
2482	setStatus(iSequence, atLowerBound);
2483	solution[iSequence] = lower[iSequence];
2484	changeObj += movement * cost[iSequence];
2485	//numberFake_--;
2486	//setFakeBound(iSequence,noFake);
2487	}
2488	#endif
2489	}
2490	}
2491	}
2492	}
2493
2494	#ifdef CLP_DEBUG
2495	if (fullRecompute && numberFake_ && (handler_->logLevel() & 16) != 0)
2496	printf("%d fake after full update\n", numberFake_);
2497	#endif
2498	// set number of infeasibilities
2499	columnArray->setNumElements(numberInfeasibilities);
2500	numberInfeasibilities += numberRowInfeasibilities;
2501	if (fullRecompute) {
2502	// do actual flips
2503	flipBounds(rowArray, columnArray);
2504	}
2505	objectiveChange += changeObj;
2506	return numberInfeasibilities;
2507	}
2508	void
2509	ClpSimplexDual::updateDualsInValuesPass(CoinIndexedVector * rowArray,
2510	CoinIndexedVector * columnArray,
2511	double theta)
2512	{
2513
2514	// use a tighter tolerance except for all being okay
2515	double tolerance = dualTolerance_;
2516
2517	// Coding is very similar but we can save a bit by splitting
2518	// Do rows
2519	{
2520	int i;
2521	double * reducedCost = djRegion(0);
2522	double * work;
2523	int number;
2524	int * which;
2525	work = rowArray->denseVector();
2526	number = rowArray->getNumElements();
2527	which = rowArray->getIndices();
2528	for (i = 0; i < number; i++) {
2529	int iSequence = which[i];
2530	double alphaI = work[i];
2531	double value = reducedCost[iSequence] - theta * alphaI;
2532	work[i] = 0.0;
2533	reducedCost[iSequence] = value;
2534
2535	Status status = getStatus(iSequence + numberColumns_);
2536	// more likely to be at upper bound ?
2537	if (status == atUpperBound) {
2538
2539	if (value > tolerance)
2540	reducedCost[iSequence] = 0.0;
2541	} else if (status == atLowerBound) {
2542
2543	if (value < -tolerance) {
2544	reducedCost[iSequence] = 0.0;
2545	}
2546	}
2547	}
2548	}
2549	rowArray->setNumElements(0);
2550
2551	// Do columns
2552	{
2553	int i;
2554	double * reducedCost = djRegion(1);
2555	double * work;
2556	int number;
2557	int * which;
2558	work = columnArray->denseVector();
2559	number = columnArray->getNumElements();
2560	which = columnArray->getIndices();
2561
2562	for (i = 0; i < number; i++) {
2563	int iSequence = which[i];
2564	double alphaI = work[i];
2565	double value = reducedCost[iSequence] - theta * alphaI;
2566	work[i] = 0.0;
2567	reducedCost[iSequence] = value;
2568
2569	Status status = getStatus(iSequence);
2570	if (status == atLowerBound) {
2571	if (value < -tolerance)
2572	reducedCost[iSequence] = 0.0;
2573	} else if (status == atUpperBound) {
2574	if (value > tolerance)
2575	reducedCost[iSequence] = 0.0;
2576	}
2577	}
2578	}
2579	columnArray->setNumElements(0);
2580	}
2581	/*
2582	Chooses dual pivot row
2583	Would be faster with separate region to scan
2584	and will have this (with square of infeasibility) when steepest
2585	For easy problems we can just choose one of the first rows we look at
2586	*/
2587	void
2588	ClpSimplexDual::dualRow(int alreadyChosen)
2589	{
2590	// get pivot row using whichever method it is
2591	int chosenRow = -1;
2592	#ifdef FORCE_FOLLOW
2593	bool forceThis = false;
2594	if (!fpFollow && strlen(forceFile)) {
2595	fpFollow = fopen(forceFile, "r");
2596	assert (fpFollow);
2597	}
2598	if (fpFollow) {
2599	if (numberIterations_ <= force_iteration) {
2600	// read to next Clp0102
2601	char temp[300];
2602	while (fgets(temp, 250, fpFollow)) {
2603	if (strncmp(temp, "Clp0102", 7))
2604	continue;
2605	char cin, cout;
2606	sscanf(temp + 9, "%d%f%s%c%c%d%s%*c%c%d",
2607	&force_iteration, &cin, &force_in, &cout, &force_out);
2608	if (cin == 'R')
2609	force_in += numberColumns_;
2610	if (cout == 'R')
2611	force_out += numberColumns_;
2612	forceThis = true;
2613	assert (numberIterations_ == force_iteration - 1);
2614	printf("Iteration %d will force %d out and %d in\n",
2615	force_iteration, force_out, force_in);
2616	alreadyChosen = force_out;
2617	break;
2618	}
2619	} else {
2620	// use old
2621	forceThis = true;
2622	}
2623	if (!forceThis) {
2624	fclose(fpFollow);
2625	fpFollow = NULL;
2626	forceFile = "";
2627	}
2628	}
2629	#endif
2630	double freeAlpha = 0.0;
2631	if (alreadyChosen < 0) {
2632	// first see if any free variables and put them in basis
2633	int nextFree = nextSuperBasic();
2634	//nextFree=-1; //off
2635	if (nextFree >= 0) {
2636	// unpack vector and find a good pivot
2637	unpack(rowArray_[1], nextFree);
2638	factorization_->updateColumn(rowArray_[2], rowArray_[1]);
2639
2640	double * work = rowArray_[1]->denseVector();
2641	int number = rowArray_[1]->getNumElements();
2642	int * which = rowArray_[1]->getIndices();
2643	double bestFeasibleAlpha = 0.0;
2644	int bestFeasibleRow = -1;
2645	double bestInfeasibleAlpha = 0.0;
2646	int bestInfeasibleRow = -1;
2647	int i;
2648
2649	for (i = 0; i < number; i++) {
2650	int iRow = which[i];
2651	double alpha = fabs(work[iRow]);
2652	if (alpha > 1.0e-3) {
2653	int iSequence = pivotVariable_[iRow];
2654	double value = solution_[iSequence];
2655	double lower = lower_[iSequence];
2656	double upper = upper_[iSequence];
2657	double infeasibility = 0.0;
2658	if (value > upper)
2659	infeasibility = value - upper;
2660	else if (value < lower)
2661	infeasibility = lower - value;
2662	if (infeasibility * alpha > bestInfeasibleAlpha && alpha > 1.0e-1) {
2663	if (!flagged(iSequence)) {
2664	bestInfeasibleAlpha = infeasibility * alpha;
2665	bestInfeasibleRow = iRow;
2666	}
2667	}
2668	if (alpha > bestFeasibleAlpha && (lower > -1.0e20 \|\| upper < 1.0e20)) {
2669	bestFeasibleAlpha = alpha;
2670	bestFeasibleRow = iRow;
2671	}
2672	}
2673	}
2674	if (bestInfeasibleRow >= 0)
2675	chosenRow = bestInfeasibleRow;
2676	else if (bestFeasibleAlpha > 1.0e-2)
2677	chosenRow = bestFeasibleRow;
2678	if (chosenRow >= 0) {
2679	pivotRow_ = chosenRow;
2680	freeAlpha = work[chosenRow];
2681	}
2682	rowArray_[1]->clear();
2683	}
2684	} else {
2685	// in values pass
2686	chosenRow = alreadyChosen;
2687	#ifdef FORCE_FOLLOW
2688	if(forceThis) {
2689	alreadyChosen = -1;
2690	chosenRow = -1;
2691	for (int i = 0; i < numberRows_; i++) {
2692	if (pivotVariable_[i] == force_out) {
2693	chosenRow = i;
2694	break;
2695	}
2696	}
2697	assert (chosenRow >= 0);
2698	}
2699	#endif
2700	pivotRow_ = chosenRow;
2701	}
2702	if (chosenRow < 0)
2703	pivotRow_ = dualRowPivot_->pivotRow();
2704
2705	if (pivotRow_ >= 0) {
2706	sequenceOut_ = pivotVariable_[pivotRow_];
2707	valueOut_ = solution_[sequenceOut_];
2708	lowerOut_ = lower_[sequenceOut_];
2709	upperOut_ = upper_[sequenceOut_];
2710	if (alreadyChosen < 0) {
2711	// if we have problems we could try other way and hope we get a
2712	// zero pivot?
2713	if (valueOut_ > upperOut_) {
2714	directionOut_ = -1;
2715	dualOut_ = valueOut_ - upperOut_;
2716	} else if (valueOut_ < lowerOut_) {
2717	directionOut_ = 1;
2718	dualOut_ = lowerOut_ - valueOut_;
2719	} else {
2720	#if 1
2721	// odd (could be free) - it's feasible - go to nearest
2722	if (valueOut_ - lowerOut_ < upperOut_ - valueOut_) {
2723	directionOut_ = 1;
2724	dualOut_ = lowerOut_ - valueOut_;
2725	} else {
2726	directionOut_ = -1;
2727	dualOut_ = valueOut_ - upperOut_;
2728	}
2729	#else
2730	// odd (could be free) - it's feasible - improve obj
2731	printf("direction from alpha of %g is %d\n",
2732	freeAlpha, freeAlpha > 0.0 ? 1 : -1);
2733	if (valueOut_ - lowerOut_ > 1.0e20)
2734	freeAlpha = 1.0;
2735	else if(upperOut_ - valueOut_ > 1.0e20)
2736	freeAlpha = -1.0;
2737	//if (valueOut_-lowerOut_<upperOut_-valueOut_) {
2738	if (freeAlpha < 0.0) {
2739	directionOut_ = 1;
2740	dualOut_ = lowerOut_ - valueOut_;
2741	} else {
2742	directionOut_ = -1;
2743	dualOut_ = valueOut_ - upperOut_;
2744	}
2745	printf("direction taken %d - bounds %g %g %g\n",
2746	directionOut_, lowerOut_, valueOut_, upperOut_);
2747	#endif
2748	}
2749	#ifdef CLP_DEBUG
2750	assert(dualOut_ >= 0.0);
2751	#endif
2752	} else {
2753	// in values pass so just use sign of dj
2754	// We don't want to go through any barriers so set dualOut low
2755	// free variables will never be here
2756	dualOut_ = 1.0e-6;
2757	if (dj_[sequenceOut_] > 0.0) {
2758	// this will give a -1 in pivot row (as slacks are -1.0)
2759	directionOut_ = 1;
2760	} else {
2761	directionOut_ = -1;
2762	}
2763	}
2764	}
2765	return ;
2766	}
2767	// Checks if any fake bounds active - if so returns number and modifies
2768	// dualBound_ and everything.
2769	// Free variables will be left as free
2770	// Returns number of bounds changed if >=0
2771	// Returns -1 if not initialize and no effect
2772	// Fills in changeVector which can be used to see if unbounded
2773	// and cost of change vector
2774	int
2775	ClpSimplexDual::changeBounds(int initialize,
2776	CoinIndexedVector * outputArray,
2777	double & changeCost)
2778	{
2779	numberFake_ = 0;
2780	if (!initialize) {
2781	int numberInfeasibilities;
2782	double newBound;
2783	newBound = 5.0 * dualBound_;
2784	numberInfeasibilities = 0;
2785	changeCost = 0.0;
2786	// put back original bounds and then check
2787	createRim1(false);
2788	int iSequence;
2789	// bounds will get bigger - just look at ones at bounds
2790	for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2791	double lowerValue = lower_[iSequence];
2792	double upperValue = upper_[iSequence];
2793	double value = solution_[iSequence];
2794	setFakeBound(iSequence, ClpSimplexDual::noFake);
2795	switch(getStatus(iSequence)) {
2796
2797	case basic:
2798	case ClpSimplex::isFixed:
2799	break;
2800	case isFree:
2801	case superBasic:
2802	break;
2803	case atUpperBound:
2804	if (fabs(value - upperValue) > primalTolerance_)
2805	numberInfeasibilities++;
2806	break;
2807	case atLowerBound:
2808	if (fabs(value - lowerValue) > primalTolerance_)
2809	numberInfeasibilities++;
2810	break;
2811	}
2812	}
2813	// If dual infeasible then carry on
2814	if (numberInfeasibilities) {
2815	handler_->message(CLP_DUAL_CHECKB, messages_)
2816	<< newBound
2817	<< CoinMessageEol;
2818	int iSequence;
2819	for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2820	double lowerValue = lower_[iSequence];
2821	double upperValue = upper_[iSequence];
2822	double newLowerValue;
2823	double newUpperValue;
2824	Status status = getStatus(iSequence);
2825	if (status == atUpperBound \|\|
2826	status == atLowerBound) {
2827	double value = solution_[iSequence];
2828	if (value - lowerValue <= upperValue - value) {
2829	newLowerValue = CoinMax(lowerValue, value - 0.666667 * newBound);
2830	newUpperValue = CoinMin(upperValue, newLowerValue + newBound);
2831	} else {
2832	newUpperValue = CoinMin(upperValue, value + 0.666667 * newBound);
2833	newLowerValue = CoinMax(lowerValue, newUpperValue - newBound);
2834	}
2835	lower_[iSequence] = newLowerValue;
2836	upper_[iSequence] = newUpperValue;
2837	if (newLowerValue > lowerValue) {
2838	if (newUpperValue < upperValue) {
2839	setFakeBound(iSequence, ClpSimplexDual::bothFake);
2840	#ifdef CLP_INVESTIGATE
2841	abort(); // No idea what should happen here - I have never got here
2842	#endif
2843	numberFake_++;
2844	} else {
2845	setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2846	numberFake_++;
2847	}
2848	} else {
2849	if (newUpperValue < upperValue) {
2850	setFakeBound(iSequence, ClpSimplexDual::upperFake);
2851	numberFake_++;
2852	}
2853	}
2854	if (status == atUpperBound)
2855	solution_[iSequence] = newUpperValue;
2856	else
2857	solution_[iSequence] = newLowerValue;
2858	double movement = solution_[iSequence] - value;
2859	if (movement && outputArray) {
2860	if (iSequence >= numberColumns_) {
2861	outputArray->quickAdd(iSequence, -movement);
2862	changeCost += movement * cost_[iSequence];
2863	} else {
2864	matrix_->add(this, outputArray, iSequence, movement);
2865	changeCost += movement * cost_[iSequence];
2866	}
2867	}
2868	}
2869	}
2870	dualBound_ = newBound;
2871	} else {
2872	numberInfeasibilities = -1;
2873	}
2874	return numberInfeasibilities;
2875	} else if (initialize == 1 \|\| initialize == 3) {
2876	int iSequence;
2877	if (initialize == 3) {
2878	for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2879	setFakeBound(iSequence, ClpSimplexDual::noFake);
2880	}
2881	}
2882	double testBound = 0.999999 * dualBound_;
2883	for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2884	Status status = getStatus(iSequence);
2885	if (status == atUpperBound \|\|
2886	status == atLowerBound) {
2887	double lowerValue = lower_[iSequence];
2888	double upperValue = upper_[iSequence];
2889	double value = solution_[iSequence];
2890	if (lowerValue > -largeValue_ \|\| upperValue < largeValue_) {
2891	if (true \|\| lowerValue - value > -0.5 * dualBound_ \|\|
2892	upperValue - value < 0.5 * dualBound_) {
2893	if (fabs(lowerValue - value) <= fabs(upperValue - value)) {
2894	if (upperValue > lowerValue + testBound) {
2895	if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2896	numberFake_++;
2897	upper_[iSequence] = lowerValue + dualBound_;
2898	setFakeBound(iSequence, ClpSimplexDual::upperFake);
2899	}
2900	} else {
2901	if (lowerValue < upperValue - testBound) {
2902	if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2903	numberFake_++;
2904	lower_[iSequence] = upperValue - dualBound_;
2905	setFakeBound(iSequence, ClpSimplexDual::lowerFake);
2906	}
2907	}
2908	} else {
2909	if (getFakeBound(iSequence) == ClpSimplexDual::noFake)
2910	numberFake_++;
2911	lower_[iSequence] = -0.5 * dualBound_;
2912	upper_[iSequence] = 0.5 * dualBound_;
2913	setFakeBound(iSequence, ClpSimplexDual::bothFake);
2914	abort();
2915	}
2916	if (status == atUpperBound)
2917	solution_[iSequence] = upper_[iSequence];
2918	else
2919	solution_[iSequence] = lower_[iSequence];
2920	} else {
2921	// set non basic free variables to fake bounds
2922	// I don't think we should ever get here
2923	CoinAssert(!("should not be here"));
2924	lower_[iSequence] = -0.5 * dualBound_;
2925	upper_[iSequence] = 0.5 * dualBound_;
2926	setFakeBound(iSequence, ClpSimplexDual::bothFake);
2927	numberFake_++;
2928	setStatus(iSequence, atUpperBound);
2929	solution_[iSequence] = 0.5 * dualBound_;
2930	}
2931	} else if (status == basic) {
2932	// make sure not at fake bound and bounds correct
2933	setFakeBound(iSequence, ClpSimplexDual::noFake);
2934	double gap = upper_[iSequence] - lower_[iSequence];
2935	if (gap > 0.5 * dualBound_ && gap < 2.0 * dualBound_) {
2936	if (iSequence < numberColumns_) {
2937	if (columnScale_) {
2938	double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2939	// lower
2940	double value = columnLower_[iSequence];
2941	if (value > -1.0e30) {
2942	value *= multiplier;
2943	}
2944	lower_[iSequence] = value;
2945	// upper
2946	value = columnUpper_[iSequence];
2947	if (value < 1.0e30) {
2948	value *= multiplier;
2949	}
2950	upper_[iSequence] = value;
2951	} else {
2952	lower_[iSequence] = columnLower_[iSequence];;
2953	upper_[iSequence] = columnUpper_[iSequence];;
2954	}
2955	} else {
2956	int iRow = iSequence - numberColumns_;
2957	if (rowScale_) {
2958	// lower
2959	double multiplier = rhsScale_ * rowScale_[iRow];
2960	double value = rowLower_[iRow];
2961	if (value > -1.0e30) {
2962	value *= multiplier;
2963	}
2964	lower_[iSequence] = value;
2965	// upper
2966	value = rowUpper_[iRow];
2967	if (value < 1.0e30) {
2968	value *= multiplier;
2969	}
2970	upper_[iSequence] = value;
2971	} else {
2972	lower_[iSequence] = rowLower_[iRow];;
2973	upper_[iSequence] = rowUpper_[iRow];;
2974	}
2975	}
2976	}
2977	}
2978	}
2979
2980	return 1;
2981	} else {
2982	// just reset changed ones
2983	if (columnScale_) {
2984	int iSequence;
2985	for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
2986	FakeBound fakeStatus = getFakeBound(iSequence);
2987	if (fakeStatus != noFake) {
2988	if ((static_cast<int> (fakeStatus) & 1) != 0) {
2989	// lower
2990	double value = columnLower_[iSequence];
2991	if (value > -1.0e30) {
2992	double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
2993	value *= multiplier;
2994	}
2995	columnLowerWork_[iSequence] = value;
2996	}
2997	if ((static_cast<int> (fakeStatus) & 2) != 0) {
2998	// upper
2999	double value = columnUpper_[iSequence];
3000	if (value < 1.0e30) {
3001	double multiplier = rhsScale_ * inverseColumnScale_[iSequence];
3002	value *= multiplier;
3003	}
3004	columnUpperWork_[iSequence] = value;
3005	}
3006	}
3007	}
3008	for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3009	FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3010	if (fakeStatus != noFake) {
3011	if ((static_cast<int> (fakeStatus) & 1) != 0) {
3012	// lower
3013	double value = rowLower_[iSequence];
3014	if (value > -1.0e30) {
3015	double multiplier = rhsScale_ * rowScale_[iSequence];
3016	value *= multiplier;
3017	}
3018	rowLowerWork_[iSequence] = value;
3019	}
3020	if ((static_cast<int> (fakeStatus) & 2) != 0) {
3021	// upper
3022	double value = rowUpper_[iSequence];
3023	if (value < 1.0e30) {
3024	double multiplier = rhsScale_ * rowScale_[iSequence];
3025	value *= multiplier;
3026	}
3027	rowUpperWork_[iSequence] = value;
3028	}
3029	}
3030	}
3031	} else {
3032	int iSequence;
3033	for (iSequence = 0; iSequence < numberColumns_; iSequence++) {
3034	FakeBound fakeStatus = getFakeBound(iSequence);
3035	if ((static_cast<int> (fakeStatus) & 1) != 0) {
3036	// lower
3037	columnLowerWork_[iSequence] = columnLower_[iSequence];
3038	}
3039	if ((static_cast<int> (fakeStatus) & 2) != 0) {
3040	// upper
3041	columnUpperWork_[iSequence] = columnUpper_[iSequence];
3042	}
3043	}
3044	for (iSequence = 0; iSequence < numberRows_; iSequence++) {
3045	FakeBound fakeStatus = getFakeBound(iSequence + numberColumns_);
3046	if ((static_cast<int> (fakeStatus) & 1) != 0) {
3047	// lower
3048	rowLowerWork_[iSequence] = rowLower_[iSequence];
3049	}
3050	if ((static_cast<int> (fakeStatus) & 2) != 0) {
3051	// upper
3052	rowUpperWork_[iSequence] = rowUpper_[iSequence];
3053	}
3054	}
3055	}
3056	return 0;
3057	}
3058	}
3059	int
3060	ClpSimplexDual::dualColumn0(const CoinIndexedVector * rowArray,
3061	const CoinIndexedVector * columnArray,
3062	CoinIndexedVector * spareArray,
3063	double acceptablePivot,
3064	double & upperReturn, double &bestReturn, double & badFree)
3065	{
3066	// do first pass to get possibles
3067	double * spare = spareArray->denseVector();
3068	int * index = spareArray->getIndices();
3069	const double * work;
3070	int number;
3071	const int * which;
3072	const double * reducedCost;
3073	// We can also see if infeasible or pivoting on free
3074	double tentativeTheta = 1.0e15;
3075	double upperTheta = 1.0e31;
3076	double freePivot = acceptablePivot;
3077	double bestPossible = 0.0;
3078	int numberRemaining = 0;
3079	int i;
3080	badFree = 0.0;
3081	if ((moreSpecialOptions_ & 8) != 0) {
3082	// No free or super basic
3083	double multiplier[] = { -1.0, 1.0};
3084	double dualT = - dualTolerance_;
3085	for (int iSection = 0; iSection < 2; iSection++) {
3086
3087	int addSequence;
3088	unsigned char * statusArray;
3089	if (!iSection) {
3090	work = rowArray->denseVector();
3091	number = rowArray->getNumElements();
3092	which = rowArray->getIndices();
3093	reducedCost = rowReducedCost_;
3094	addSequence = numberColumns_;
3095	statusArray = status_ + numberColumns_;
3096	} else {
3097	work = columnArray->denseVector();
3098	number = columnArray->getNumElements();
3099	which = columnArray->getIndices();
3100	reducedCost = reducedCostWork_;
3101	addSequence = 0;
3102	statusArray = status_;
3103	}
3104
3105	for (i = 0; i < number; i++) {
3106	int iSequence = which[i];
3107	double alpha;
3108	double oldValue;
3109	double value;
3110
3111	assert (getStatus(iSequence + addSequence) != isFree
3112	&& getStatus(iSequence + addSequence) != superBasic);
3113	int iStatus = (statusArray[iSequence] & 3) - 1;
3114	if (iStatus) {
3115	double mult = multiplier[iStatus-1];
3116	alpha = work[i] * mult;
3117	if (alpha > 0.0) {
3118	oldValue = reducedCost[iSequence] * mult;
3119	value = oldValue - tentativeTheta * alpha;
3120	if (value < dualT) {
3121	bestPossible = CoinMax(bestPossible, alpha);
3122	value = oldValue - upperTheta * alpha;
3123	if (value < dualT && alpha >= acceptablePivot) {
3124	upperTheta = (oldValue - dualT) / alpha;
3125	//tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3126	}
3127	// add to list
3128	spare[numberRemaining] = alpha * mult;
3129	index[numberRemaining++] = iSequence + addSequence;
3130	}
3131	}
3132	}
3133	}
3134	}
3135	} else {
3136	// some free or super basic
3137	for (int iSection = 0; iSection < 2; iSection++) {
3138
3139	int addSequence;
3140
3141	if (!iSection) {
3142	work = rowArray->denseVector();
3143	number = rowArray->getNumElements();
3144	which = rowArray->getIndices();
3145	reducedCost = rowReducedCost_;
3146	addSequence = numberColumns_;
3147	} else {
3148	work = columnArray->denseVector();
3149	number = columnArray->getNumElements();
3150	which = columnArray->getIndices();
3151	reducedCost = reducedCostWork_;
3152	addSequence = 0;
3153	}
3154
3155	for (i = 0; i < number; i++) {
3156	int iSequence = which[i];
3157	double alpha;
3158	double oldValue;
3159	double value;
3160	bool keep;
3161
3162	switch(getStatus(iSequence + addSequence)) {
3163
3164	case basic:
3165	case ClpSimplex::isFixed:
3166	break;
3167	case isFree:
3168	case superBasic:
3169	alpha = work[i];
3170	bestPossible = CoinMax(bestPossible, fabs(alpha));
3171	oldValue = reducedCost[iSequence];
3172	// If free has to be very large - should come in via dualRow
3173	//if (getStatus(iSequence+addSequence)==isFree&&fabs(alpha)<1.0e-3)
3174	//break;
3175	if (oldValue > dualTolerance_) {
3176	keep = true;
3177	} else if (oldValue < -dualTolerance_) {
3178	keep = true;
3179	} else {
3180	if (fabs(alpha) > CoinMax(10.0 * acceptablePivot, 1.0e-5)) {
3181	keep = true;
3182	} else {
3183	keep = false;
3184	badFree = CoinMax(badFree, fabs(alpha));
3185	}
3186	}
3187	if (keep) {
3188	// free - choose largest
3189	if (fabs(alpha) > freePivot) {
3190	freePivot = fabs(alpha);
3191	sequenceIn_ = iSequence + addSequence;
3192	theta_ = oldValue / alpha;
3193	alpha_ = alpha;
3194	}
3195	}
3196	break;
3197	case atUpperBound:
3198	alpha = work[i];
3199	oldValue = reducedCost[iSequence];
3200	value = oldValue - tentativeTheta * alpha;
3201	//assert (oldValue<=dualTolerance_*1.0001);
3202	if (value > dualTolerance_) {
3203	bestPossible = CoinMax(bestPossible, -alpha);
3204	value = oldValue - upperTheta * alpha;
3205	if (value > dualTolerance_ && -alpha >= acceptablePivot) {
3206	upperTheta = (oldValue - dualTolerance_) / alpha;
3207	//tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3208	}
3209	// add to list
3210	spare[numberRemaining] = alpha;
3211	index[numberRemaining++] = iSequence + addSequence;
3212	}
3213	break;
3214	case atLowerBound:
3215	alpha = work[i];
3216	oldValue = reducedCost[iSequence];
3217	value = oldValue - tentativeTheta * alpha;
3218	//assert (oldValue>=-dualTolerance_*1.0001);
3219	if (value < -dualTolerance_) {
3220	bestPossible = CoinMax(bestPossible, alpha);
3221	value = oldValue - upperTheta * alpha;
3222	if (value < -dualTolerance_ && alpha >= acceptablePivot) {
3223	upperTheta = (oldValue + dualTolerance_) / alpha;
3224	//tentativeTheta = CoinMin(2.0*upperTheta,tentativeTheta);
3225	}
3226	// add to list
3227	spare[numberRemaining] = alpha;
3228	index[numberRemaining++] = iSequence + addSequence;
3229	}
3230	break;
3231	}
3232	}
3233	}
3234	}
3235	upperReturn = upperTheta;
3236	bestReturn = bestPossible;
3237	return numberRemaining;
3238	}
3239	/*
3240	Row array has row part of pivot row (as duals so sign may be switched)
3241	Column array has column part.
3242	This chooses pivot column.
3243	Spare array will be needed when we start getting clever.
3244	We will check for basic so spare array will never overflow.
3245	If necessary will modify costs
3246	*/
3247	double
3248	ClpSimplexDual::dualColumn(CoinIndexedVector * rowArray,
3249	CoinIndexedVector * columnArray,
3250	CoinIndexedVector * spareArray,
3251	CoinIndexedVector * spareArray2,
3252	double acceptablePivot,
3253	CoinBigIndex * /dubiousWeights/)
3254	{
3255	int numberPossiblySwapped = 0;
3256	int numberRemaining = 0;
3257
3258	double totalThru = 0.0; // for when variables flip
3259	//double saveAcceptable=acceptablePivot;
3260	//acceptablePivot=1.0e-9;
3261
3262	double bestEverPivot = acceptablePivot;
3263	int lastSequence = -1;
3264	double lastPivot = 0.0;
3265	double upperTheta;
3266	double newTolerance = dualTolerance_;
3267	//newTolerance = dualTolerance_+1.0e-6*dblParam_[ClpDualTolerance];
3268	// will we need to increase tolerance
3269	bool thisIncrease = false;
3270	// If we think we need to modify costs (not if something from broad sweep)
3271	bool modifyCosts = false;
3272	// Increase in objective due to swapping bounds (may be negative)
3273	double increaseInObjective = 0.0;
3274
3275	// use spareArrays to put ones looked at in
3276	// we are going to flip flop between
3277	int iFlip = 0;
3278	// Possible list of pivots
3279	int interesting[2];
3280	// where possible swapped ones are
3281	int swapped[2];
3282	// for zeroing out arrays after
3283	int marker[2][2];
3284	// pivot elements
3285	double * array[2], * spare, * spare2;
3286	// indices
3287	int * indices[2], * index, * index2;
3288	spareArray2->clear();
3289	array[0] = spareArray->denseVector();
3290	indices[0] = spareArray->getIndices();
3291	spare = array[0];
3292	index = indices[0];
3293	array[1] = spareArray2->denseVector();
3294	indices[1] = spareArray2->getIndices();
3295	int i;
3296
3297	// initialize lists
3298	for (i = 0; i < 2; i++) {
3299	interesting[i] = 0;
3300	swapped[i] = numberColumns_;
3301	marker[i][0] = 0;
3302	marker[i][1] = numberColumns_;
3303	}
3304	/*
3305	First we get a list of possible pivots. We can also see if the
3306	problem looks infeasible or whether we want to pivot in free variable.
3307	This may make objective go backwards but can only happen a finite
3308	number of times and I do want free variables basic.
3309
3310	Then we flip back and forth. At the start of each iteration
3311	interesting[iFlip] should have possible candidates and swapped[iFlip]
3312	will have pivots if we decide to take a previous pivot.
3313	At end of each iteration interesting[1-iFlip] should have
3314	candidates if we go through this theta and swapped[1-iFlip]
3315	pivots if we don't go through.
3316
3317	At first we increase theta and see what happens. We start
3318	theta at a reasonable guess. If in right area then we do bit by bit.
3319
3320	*/
3321
3322	// do first pass to get possibles
3323	upperTheta = 1.0e31;
3324	double bestPossible = 0.0;
3325	double badFree = 0.0;
3326	alpha_ = 0.0;
3327	if (spareIntArray_[0] >= 0) {
3328	numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3329	acceptablePivot, upperTheta, bestPossible, badFree);
3330	} else {
3331	// already done
3332	numberRemaining = spareArray->getNumElements();
3333	spareArray->setNumElements(0);
3334	upperTheta = spareDoubleArray_[0];
3335	bestPossible = spareDoubleArray_[1];
3336	if (spareIntArray_[0] == -1) {
3337	theta_ = spareDoubleArray_[2];
3338	alpha_ = spareDoubleArray_[3];
3339	sequenceIn_ = spareIntArray_[1];
3340	} else {
3341	#if 0
3342	int n = numberRemaining;
3343	double u = upperTheta;
3344	double b = bestPossible;
3345	upperTheta = 1.0e31;
3346	bestPossible = 0.0;
3347	numberRemaining = dualColumn0(rowArray, columnArray, spareArray,
3348	acceptablePivot, upperTheta, bestPossible, badFree);
3349	assert (n == numberRemaining);
3350	assert (fabs(b - bestPossible) < 1.0e-7);
3351	assert (fabs(u - upperTheta) < 1.0e-7);
3352	#endif
3353	}
3354	}
3355	// switch off
3356	spareIntArray_[0] = 0;
3357	// We can also see if infeasible or pivoting on free
3358	double tentativeTheta = 1.0e25;
3359	interesting[0] = numberRemaining;
3360	marker[0][0] = numberRemaining;
3361
3362	if (!numberRemaining && sequenceIn_ < 0)
3363	return 0.0; // Looks infeasible
3364
3365	// If sum of bad small pivots too much
3366	#define MORE_CAREFUL
3367	#ifdef MORE_CAREFUL
3368	bool badSumPivots = false;
3369	#endif
3370	if (sequenceIn_ >= 0) {
3371	// free variable - always choose
3372	} else {
3373
3374	theta_ = 1.0e50;
3375	// now flip flop between spare arrays until reasonable theta
3376	tentativeTheta = CoinMax(10.0 * upperTheta, 1.0e-7);
3377
3378	// loops increasing tentative theta until can't go through
3379
3380	while (tentativeTheta < 1.0e22) {
3381	double thruThis = 0.0;
3382
3383	double bestPivot = acceptablePivot;
3384	int bestSequence = -1;
3385
3386	numberPossiblySwapped = numberColumns_;
3387	numberRemaining = 0;
3388
3389	upperTheta = 1.0e50;
3390
3391	spare = array[iFlip];
3392	index = indices[iFlip];
3393	spare2 = array[1-iFlip];
3394	index2 = indices[1-iFlip];
3395
3396	// try 3 different ways
3397	// 1 bias increase by ones with slightly wrong djs
3398	// 2 bias by all
3399	// 3 bias by all - tolerance
3400	#define TRYBIAS 3
3401
3402
3403	double increaseInThis = 0.0; //objective increase in this loop
3404
3405	for (i = 0; i < interesting[iFlip]; i++) {
3406	int iSequence = index[i];
3407	double alpha = spare[i];
3408	double oldValue = dj_[iSequence];
3409	double value = oldValue - tentativeTheta * alpha;
3410
3411	if (alpha < 0.0) {
3412	//at upper bound
3413	if (value > newTolerance) {
3414	double range = upper_[iSequence] - lower_[iSequence];
3415	thruThis -= range * alpha;
3416	#if TRYBIAS==1
3417	if (oldValue > 0.0)
3418	increaseInThis -= oldValue * range;
3419	#elif TRYBIAS==2
3420	increaseInThis -= oldValue * range;
3421	#else
3422	increaseInThis -= (oldValue + dualTolerance_) * range;
3423	#endif
3424	// goes on swapped list (also means candidates if too many)
3425	spare2[--numberPossiblySwapped] = alpha;
3426	index2[numberPossiblySwapped] = iSequence;
3427	if (fabs(alpha) > bestPivot) {
3428	bestPivot = fabs(alpha);
3429	bestSequence = numberPossiblySwapped;
3430	}
3431	} else {
3432	value = oldValue - upperTheta * alpha;
3433	if (value > newTolerance && -alpha >= acceptablePivot)
3434	upperTheta = (oldValue - newTolerance) / alpha;
3435	spare2[numberRemaining] = alpha;
3436	index2[numberRemaining++] = iSequence;
3437	}
3438	} else {
3439	// at lower bound
3440	if (value < -newTolerance) {
3441	double range = upper_[iSequence] - lower_[iSequence];
3442	thruThis += range * alpha;
3443	//?? is this correct - and should we look at good ones
3444	#if TRYBIAS==1
3445	if (oldValue < 0.0)
3446	increaseInThis += oldValue * range;
3447	#elif TRYBIAS==2
3448	increaseInThis += oldValue * range;
3449	#else
3450	increaseInThis += (oldValue - dualTolerance_) * range;
3451	#endif
3452	// goes on swapped list (also means candidates if too many)
3453	spare2[--numberPossiblySwapped] = alpha;
3454	index2[numberPossiblySwapped] = iSequence;
3455	if (fabs(alpha) > bestPivot) {
3456	bestPivot = fabs(alpha);
3457	bestSequence = numberPossiblySwapped;
3458	}
3459	} else {
3460	value = oldValue - upperTheta * alpha;
3461	if (value < -newTolerance && alpha >= acceptablePivot)
3462	upperTheta = (oldValue + newTolerance) / alpha;
3463	spare2[numberRemaining] = alpha;
3464	index2[numberRemaining++] = iSequence;
3465	}
3466	}
3467	}
3468	swapped[1-iFlip] = numberPossiblySwapped;
3469	interesting[1-iFlip] = numberRemaining;
3470	marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3471	marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3472
3473	if (totalThru + thruThis >= fabs(dualOut_) \|\|
3474	increaseInObjective + increaseInThis < 0.0) {
3475	// We should be pivoting in this batch
3476	// so compress down to this lot
3477	numberRemaining = 0;
3478	for (i = numberColumns_ - 1; i >= swapped[1-iFlip]; i--) {
3479	spare[numberRemaining] = spare2[i];
3480	index[numberRemaining++] = index2[i];
3481	}
3482	interesting[iFlip] = numberRemaining;
3483	int iTry;
3484	#define MAXTRY 100
3485	// first get ratio with tolerance
3486	for (iTry = 0; iTry < MAXTRY; iTry++) {
3487
3488	upperTheta = 1.0e50;
3489	numberPossiblySwapped = numberColumns_;
3490	numberRemaining = 0;
3491
3492	increaseInThis = 0.0; //objective increase in this loop
3493
3494	thruThis = 0.0;
3495
3496	spare = array[iFlip];
3497	index = indices[iFlip];
3498	spare2 = array[1-iFlip];
3499	index2 = indices[1-iFlip];
3500	for (i = 0; i < interesting[iFlip]; i++) {
3501	int iSequence = index[i];
3502	double alpha = spare[i];
3503	double oldValue = dj_[iSequence];
3504	double value = oldValue - upperTheta * alpha;
3505
3506	if (alpha < 0.0) {
3507	//at upper bound
3508	if (value > newTolerance) {
3509	if (-alpha >= acceptablePivot) {
3510	upperTheta = (oldValue - newTolerance) / alpha;
3511	}
3512	}
3513	} else {
3514	// at lower bound
3515	if (value < -newTolerance) {
3516	if (alpha >= acceptablePivot) {
3517	upperTheta = (oldValue + newTolerance) / alpha;
3518	}
3519	}
3520	}
3521	}
3522	bestPivot = acceptablePivot;
3523	sequenceIn_ = -1;
3524	#ifdef DUBIOUS_WEIGHTS
3525	double bestWeight = COIN_DBL_MAX;
3526	#endif
3527	double largestPivot = acceptablePivot;
3528	// now choose largest and sum all ones which will go through
3529	//printf("XX it %d number %d\n",numberIterations_,interesting[iFlip]);
3530	// Sum of bad small pivots
3531	#ifdef MORE_CAREFUL
3532	double sumBadPivots = 0.0;
3533	badSumPivots = false;
3534	#endif
3535	// Make sure upperTheta will work (-O2 and above gives problems)
3536	upperTheta *= 1.0000000001;
3537	for (i = 0; i < interesting[iFlip]; i++) {
3538	int iSequence = index[i];
3539	double alpha = spare[i];
3540	double value = dj_[iSequence] - upperTheta * alpha;
3541	double badDj = 0.0;
3542
3543	bool addToSwapped = false;
3544
3545	if (alpha < 0.0) {
3546	//at upper bound
3547	if (value >= 0.0) {
3548	addToSwapped = true;
3549	#if TRYBIAS==1
3550	badDj = -CoinMax(dj_[iSequence], 0.0);
3551	#elif TRYBIAS==2
3552	badDj = -dj_[iSequence];
3553	#else
3554	badDj = -dj_[iSequence] - dualTolerance_;
3555	#endif
3556	}
3557	} else {
3558	// at lower bound
3559	if (value <= 0.0) {
3560	addToSwapped = true;
3561	#if TRYBIAS==1
3562	badDj = CoinMin(dj_[iSequence], 0.0);
3563	#elif TRYBIAS==2
3564	badDj = dj_[iSequence];
3565	#else
3566	badDj = dj_[iSequence] - dualTolerance_;
3567	#endif
3568	}
3569	}
3570	if (!addToSwapped) {
3571	// add to list of remaining
3572	spare2[numberRemaining] = alpha;
3573	index2[numberRemaining++] = iSequence;
3574	} else {
3575	// add to list of swapped
3576	spare2[--numberPossiblySwapped] = alpha;
3577	index2[numberPossiblySwapped] = iSequence;
3578	// select if largest pivot
3579	bool take = false;
3580	double absAlpha = fabs(alpha);
3581	#ifdef DUBIOUS_WEIGHTS
3582	// User could do anything to break ties here
3583	double weight;
3584	if (dubiousWeights)
3585	weight = dubiousWeights[iSequence];
3586	else
3587	weight = 1.0;
3588	weight += randomNumberGenerator_.randomDouble() * 1.0e-2;
3589	if (absAlpha > 2.0 * bestPivot) {
3590	take = true;
3591	} else if (absAlpha > largestPivot) {
3592	// could multiply absAlpha and weight
3593	if (weight * bestPivot < bestWeight * absAlpha)
3594	take = true;
3595	}
3596	#else
3597	if (absAlpha > bestPivot)
3598	take = true;
3599	#endif
3600	#ifdef MORE_CAREFUL
3601	if (absAlpha < acceptablePivot && upperTheta < 1.0e20) {
3602	if (alpha < 0.0) {
3603	//at upper bound
3604	if (value > dualTolerance_) {
3605	double gap = upper_[iSequence] - lower_[iSequence];
3606	if (gap < 1.0e20)
3607	sumBadPivots += value * gap;
3608	else
3609	sumBadPivots += 1.0e20;
3610	//printf("bad %d alpha %g dj at upper %g\n",
3611	// iSequence,alpha,value);
3612	}
3613	} else {
3614	//at lower bound
3615	if (value < -dualTolerance_) {
3616	double gap = upper_[iSequence] - lower_[iSequence];
3617	if (gap < 1.0e20)
3618	sumBadPivots -= value * gap;
3619	else
3620	sumBadPivots += 1.0e20;
3621	//printf("bad %d alpha %g dj at lower %g\n",
3622	// iSequence,alpha,value);
3623	}
3624	}
3625	}
3626	#endif
3627	#ifdef FORCE_FOLLOW
3628	if (iSequence == force_in) {
3629	printf("taking %d - alpha %g best %g\n", force_in, absAlpha, largestPivot);
3630	take = true;
3631	}
3632	#endif
3633	if (take) {
3634	sequenceIn_ = numberPossiblySwapped;
3635	bestPivot = absAlpha;
3636	theta_ = dj_[iSequence] / alpha;
3637	largestPivot = CoinMax(largestPivot, 0.5 * bestPivot);
3638	#ifdef DUBIOUS_WEIGHTS
3639	bestWeight = weight;
3640	#endif
3641	//printf(" taken seq %d alpha %g weight %d\n",
3642	// iSequence,absAlpha,dubiousWeights[iSequence]);
3643	} else {
3644	//printf(" not taken seq %d alpha %g weight %d\n",
3645	// iSequence,absAlpha,dubiousWeights[iSequence]);
3646	}
3647	double range = upper_[iSequence] - lower_[iSequence];
3648	thruThis += range * fabs(alpha);
3649	increaseInThis += badDj * range;
3650	}
3651	}
3652	marker[1-iFlip][0] = CoinMax(marker[1-iFlip][0], numberRemaining);
3653	marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], numberPossiblySwapped);
3654	#ifdef MORE_CAREFUL
3655	// If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3656	if (sumBadPivots > 1.0e4) {
3657	if (handler_->logLevel() > 1)
3658	printf("maybe forcing re-factorization - sum %g %d pivots\n", sumBadPivots,
3659	factorization_->pivots());
3660	if(factorization_->pivots() > 3) {
3661	badSumPivots = true;
3662	break;
3663	}
3664	}
3665	#endif
3666	swapped[1-iFlip] = numberPossiblySwapped;
3667	interesting[1-iFlip] = numberRemaining;
3668	// If we stop now this will be increase in objective (I think)
3669	double increase = (fabs(dualOut_) - totalThru) * theta_;
3670	increase += increaseInObjective;
3671	if (theta_ < 0.0)
3672	thruThis += fabs(dualOut_); // force using this one
3673	if (increaseInObjective < 0.0 && increase < 0.0 && lastSequence >= 0) {
3674	// back
3675	// We may need to be more careful - we could do by
3676	// switch so we always do fine grained?
3677	bestPivot = 0.0;
3678	} else {
3679	// add in
3680	totalThru += thruThis;
3681	increaseInObjective += increaseInThis;
3682	}
3683	if (bestPivot < 0.1 * bestEverPivot &&
3684	bestEverPivot > 1.0e-6 &&
3685	(bestPivot < 1.0e-3 \|\| totalThru * 2.0 > fabs(dualOut_))) {
3686	// back to previous one
3687	sequenceIn_ = lastSequence;
3688	// swap regions
3689	iFlip = 1 - iFlip;
3690	break;
3691	} else if (sequenceIn_ == -1 && upperTheta > largeValue_) {
3692	if (lastPivot > acceptablePivot) {
3693	// back to previous one
3694	sequenceIn_ = lastSequence;
3695	// swap regions
3696	iFlip = 1 - iFlip;
3697	} else {
3698	// can only get here if all pivots too small
3699	}
3700	break;
3701	} else if (totalThru >= fabs(dualOut_)) {
3702	modifyCosts = true; // fine grain - we can modify costs
3703	break; // no point trying another loop
3704	} else {
3705	lastSequence = sequenceIn_;
3706	if (bestPivot > bestEverPivot)
3707	bestEverPivot = bestPivot;
3708	iFlip = 1 - iFlip;
3709	modifyCosts = true; // fine grain - we can modify costs
3710	}
3711	}
3712	if (iTry == MAXTRY)
3713	iFlip = 1 - iFlip; // flip back
3714	break;
3715	} else {
3716	// skip this lot
3717	if (bestPivot > 1.0e-3 \|\| bestPivot > bestEverPivot) {
3718	bestEverPivot = bestPivot;
3719	lastSequence = bestSequence;
3720	} else {
3721	// keep old swapped
3722	CoinMemcpyN(array[iFlip] + swapped[iFlip],
3723	numberColumns_ - swapped[iFlip], array[1-iFlip] + swapped[iFlip]);
3724	CoinMemcpyN(indices[iFlip] + swapped[iFlip],
3725	numberColumns_ - swapped[iFlip], indices[1-iFlip] + swapped[iFlip]);
3726	marker[1-iFlip][1] = CoinMin(marker[1-iFlip][1], swapped[iFlip]);
3727	swapped[1-iFlip] = swapped[iFlip];
3728	}
3729	increaseInObjective += increaseInThis;
3730	iFlip = 1 - iFlip; // swap regions
3731	tentativeTheta = 2.0 * upperTheta;
3732	totalThru += thruThis;
3733	}
3734	}
3735
3736	// can get here without sequenceIn_ set but with lastSequence
3737	if (sequenceIn_ < 0 && lastSequence >= 0) {
3738	// back to previous one
3739	sequenceIn_ = lastSequence;
3740	// swap regions
3741	iFlip = 1 - iFlip;
3742	}
3743
3744	#define MINIMUMTHETA 1.0e-18
3745	// Movement should be minimum for anti-degeneracy - unless
3746	// fixed variable out
3747	double minimumTheta;
3748	if (upperOut_ > lowerOut_)
3749	minimumTheta = MINIMUMTHETA;
3750	else
3751	minimumTheta = 0.0;
3752	if (sequenceIn_ >= 0) {
3753	// at this stage sequenceIn_ is just pointer into index array
3754	// flip just so we can use iFlip
3755	iFlip = 1 - iFlip;
3756	spare = array[iFlip];
3757	index = indices[iFlip];
3758	double oldValue;
3759	alpha_ = spare[sequenceIn_];
3760	sequenceIn_ = indices[iFlip][sequenceIn_];
3761	oldValue = dj_[sequenceIn_];
3762	theta_ = CoinMax(oldValue / alpha_, 0.0);
3763	if (theta_ < minimumTheta && fabs(alpha_) < 1.0e5 && 1) {
3764	// can't pivot to zero
3765	#if 0
3766	if (oldValue - minimumTheta*alpha_ >= -dualTolerance_) {
3767	theta_ = minimumTheta;
3768	} else if (oldValue - minimumTheta*alpha_ >= -newTolerance) {
3769	theta_ = minimumTheta;
3770	thisIncrease = true;
3771	} else {
3772	theta_ = CoinMax((oldValue + newTolerance) / alpha_, 0.0);
3773	thisIncrease = true;
3774	}
3775	#else
3776	theta_ = minimumTheta;
3777	#endif
3778	}
3779	// may need to adjust costs so all dual feasible AND pivoted is exactly 0
3780	//int costOffset = numberRows_+numberColumns_;
3781	if (modifyCosts && !badSumPivots) {
3782	int i;
3783	for (i = numberColumns_ - 1; i >= swapped[iFlip]; i--) {
3784	int iSequence = index[i];
3785	double alpha = spare[i];
3786	double value = dj_[iSequence] - theta_ * alpha;
3787
3788	// can't be free here
3789
3790	if (alpha < 0.0) {
3791	//at upper bound
3792	if (value > dualTolerance_) {
3793	thisIncrease = true;
3794	#if CLP_CAN_HAVE_ZERO_OBJ<2
3795	#define MODIFYCOST 2
3796	#endif
3797	#if MODIFYCOST
3798	// modify cost to hit new tolerance
3799	double modification = alpha * theta_ - dj_[iSequence]
3800	+ newTolerance;
3801	if ((specialOptions_&(2048 + 4096 + 16384)) != 0) {
3802	if ((specialOptions_ & 16384) != 0) {
3803	if (fabs(modification) < 1.0e-8)
3804	modification = 0.0;
3805	} else if ((specialOptions_ & 2048) != 0) {
3806	if (fabs(modification) < 1.0e-10)
3807	modification = 0.0;
3808	} else {
3809	if (fabs(modification) < 1.0e-12)
3810	modification = 0.0;
3811	}
3812	}
3813	dj_[iSequence] += modification;
3814	cost_[iSequence] += modification;
3815	if (modification)
3816	numberChanged_ ++; // Say changed costs
3817	//cost_[iSequence+costOffset] += modification; // save change
3818	#endif
3819	}
3820	} else {
3821	// at lower bound
3822	if (-value > dualTolerance_) {
3823	thisIncrease = true;
3824	#if MODIFYCOST
3825	// modify cost to hit new tolerance
3826	double modification = alpha * theta_ - dj_[iSequence]
3827	- newTolerance;
3828	//modification = CoinMax(modification,-dualTolerance_);
3829	//assert (fabs(modification)<1.0e-7);
3830	if ((specialOptions_&(2048 + 4096)) != 0) {
3831	if ((specialOptions_ & 2048) != 0) {
3832	if (fabs(modification) < 1.0e-10)
3833	modification = 0.0;
3834	} else {
3835	if (fabs(modification) < 1.0e-12)
3836	modification = 0.0;
3837	}
3838	}
3839	dj_[iSequence] += modification;
3840	cost_[iSequence] += modification;
3841	if (modification)
3842	numberChanged_ ++; // Say changed costs
3843	//cost_[iSequence+costOffset] += modification; // save change
3844	#endif
3845	}
3846	}
3847	}
3848	}
3849	}
3850	}
3851
3852	#ifdef MORE_CAREFUL
3853	// If we have done pivots and things look bad set alpha_ 0.0 to force factorization
3854	if ((badSumPivots \|\|
3855	fabs(theta_ * badFree) > 10.0 * dualTolerance_) && factorization_->pivots()) {
3856	if (handler_->logLevel() > 1)
3857	printf("forcing re-factorization\n");
3858	sequenceIn_ = -1;
3859	}
3860	#endif
3861	if (sequenceIn_ >= 0) {
3862	lowerIn_ = lower_[sequenceIn_];
3863	upperIn_ = upper_[sequenceIn_];
3864	valueIn_ = solution_[sequenceIn_];
3865	dualIn_ = dj_[sequenceIn_];
3866
3867	if (numberTimesOptimal_) {
3868	// can we adjust cost back closer to original
3869	//*** add coding
3870	}
3871	#if MODIFYCOST>1
3872	// modify cost to hit zero exactly
3873	// so (dualIn_+modification)==theta_*alpha_
3874	double modification = theta_ * alpha_ - dualIn_;
3875	// But should not move objective too much ??
3876	#define DONT_MOVE_OBJECTIVE
3877	#ifdef DONT_MOVE_OBJECTIVE
3878	double moveObjective = fabs(modification * solution_[sequenceIn_]);
3879	double smallMove = CoinMax(fabs(objectiveValue_), 1.0e-3);
3880	if (moveObjective > smallMove) {
3881	if (handler_->logLevel() > 1)
3882	printf("would move objective by %g - original mod %g sol value %g\n", moveObjective,
3883	modification, solution_[sequenceIn_]);
3884	modification *= smallMove / moveObjective;
3885	}
3886	#endif
3887	if (badSumPivots)
3888	modification = 0.0;
3889	if ((specialOptions_&(2048 + 4096)) != 0) {
3890	if ((specialOptions_ & 16384) != 0) {
3891	// in fast dual
3892	if (fabs(modification) < 1.0e-7)
3893	modification = 0.0;
3894	} else if ((specialOptions_ & 2048) != 0) {
3895	if (fabs(modification) < 1.0e-10)
3896	modification = 0.0;
3897	} else {
3898	if (fabs(modification) < 1.0e-12)
3899	modification = 0.0;
3900	}
3901	}
3902	dualIn_ += modification;
3903	dj_[sequenceIn_] = dualIn_;
3904	cost_[sequenceIn_] += modification;
3905	if (modification)
3906	numberChanged_ ++; // Say changed costs
3907	//int costOffset = numberRows_+numberColumns_;
3908	//cost_[sequenceIn_+costOffset] += modification; // save change
3909	//assert (fabs(modification)<1.0e-6);
3910	#ifdef CLP_DEBUG
3911	if ((handler_->logLevel() & 32) && fabs(modification) > 1.0e-15)
3912	printf("exact %d new cost %g, change %g\n", sequenceIn_,
3913	cost_[sequenceIn_], modification);
3914	#endif
3915	#endif
3916
3917	if (alpha_ < 0.0) {
3918	// as if from upper bound
3919	directionIn_ = -1;
3920	upperIn_ = valueIn_;
3921	} else {
3922	// as if from lower bound
3923	directionIn_ = 1;
3924	lowerIn_ = valueIn_;
3925	}
3926	} else {
3927	// no pivot
3928	bestPossible = 0.0;
3929	alpha_ = 0.0;
3930	}
3931	//if (thisIncrease)
3932	//dualTolerance_+= 1.0e-6*dblParam_[ClpDualTolerance];
3933
3934	// clear arrays
3935
3936	for (i = 0; i < 2; i++) {
3937	CoinZeroN(array[i], marker[i][0]);
3938	CoinZeroN(array[i] + marker[i][1], numberColumns_ - marker[i][1]);
3939	}
3940	return bestPossible;
3941	}
3942	#ifdef CLP_ALL_ONE_FILE
3943	#undef MAXTRY
3944	#endif
3945	/* Checks if tentative optimal actually means unbounded
3946	Returns -3 if not, 2 if is unbounded */
3947	int
3948	ClpSimplexDual::checkUnbounded(CoinIndexedVector * ray,
3949	CoinIndexedVector * spare,
3950	double changeCost)
3951	{
3952	int status = 2; // say unbounded
3953	factorization_->updateColumn(spare, ray);
3954	// get reduced cost
3955	int i;
3956	int number = ray->getNumElements();
3957	int * index = ray->getIndices();
3958	double * array = ray->denseVector();
3959	for (i = 0; i < number; i++) {
3960	int iRow = index[i];
3961	int iPivot = pivotVariable_[iRow];
3962	changeCost -= cost(iPivot) * array[iRow];
3963	}
3964	double way;
3965	if (changeCost > 0.0) {
3966	//try going down
3967	way = 1.0;
3968	} else if (changeCost < 0.0) {
3969	//try going up
3970	way = -1.0;
3971	} else {
3972	#ifdef CLP_DEBUG
3973	printf("can't decide on up or down\n");
3974	#endif
3975	way = 0.0;
3976	status = -3;
3977	}
3978	double movement = 1.0e10 * way; // some largish number
3979	double zeroTolerance = 1.0e-14 * dualBound_;
3980	for (i = 0; i < number; i++) {
3981	int iRow = index[i];
3982	int iPivot = pivotVariable_[iRow];
3983	double arrayValue = array[iRow];
3984	if (fabs(arrayValue) < zeroTolerance)
3985	arrayValue = 0.0;
3986	double newValue = solution(iPivot) + movement * arrayValue;
3987	if (newValue > upper(iPivot) + primalTolerance_ \|\|
3988	newValue < lower(iPivot) - primalTolerance_)
3989	status = -3; // not unbounded
3990	}
3991	if (status == 2) {
3992	// create ray
3993	delete [] ray_;
3994	ray_ = new double [numberColumns_];
3995	CoinZeroN(ray_, numberColumns_);
3996	for (i = 0; i < number; i++) {
3997	int iRow = index[i];
3998	int iPivot = pivotVariable_[iRow];
3999	double arrayValue = array[iRow];
4000	if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
4001	ray_[iPivot] = way * array[iRow];
4002	}
4003	}
4004	ray->clear();
4005	return status;
4006	}
4007	//static int count_alpha=0;
4008	/* Checks if finished. Updates status */
4009	void
4010	ClpSimplexDual::statusOfProblemInDual(int & lastCleaned, int type,
4011	double * givenDuals, ClpDataSave & saveData,
4012	int ifValuesPass)
4013	{
4014	#ifdef CLP_INVESTIGATE_SERIAL
4015	if (z_thinks > 0 && z_thinks < 2)
4016	z_thinks += 2;
4017	#endif
4018	bool arraysNotCreated = (type==0);
4019	// If lots of iterations then adjust costs if large ones
4020	if (numberIterations_ > 4 * (numberRows_ + numberColumns_) && objectiveScale_ == 1.0) {
4021	double largest = 0.0;
4022	for (int i = 0; i < numberRows_; i++) {
4023	int iColumn = pivotVariable_[i];
4024	largest = CoinMax(largest, fabs(cost_[iColumn]));
4025	}
4026	if (largest > 1.0e6) {
4027	objectiveScale_ = 1.0e6 / largest;
4028	for (int i = 0; i < numberRows_ + numberColumns_; i++)
4029	cost_[i] *= objectiveScale_;
4030	}
4031	}
4032	int numberPivots = factorization_->pivots();
4033	double realDualInfeasibilities = 0.0;
4034	if (type == 2) {
4035	if (alphaAccuracy_ != -1.0)
4036	alphaAccuracy_ = -2.0;
4037	// trouble - restore solution
4038	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4039	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4040	numberRows_, rowActivityWork_);
4041	CoinMemcpyN(savedSolution_ ,
4042	numberColumns_, columnActivityWork_);
4043	// restore extra stuff
4044	int dummy;
4045	matrix_->generalExpanded(this, 6, dummy);
4046	forceFactorization_ = 1; // a bit drastic but ..
4047	changeMade_++; // say something changed
4048	// get correct bounds on all variables
4049	resetFakeBounds(0);
4050	}
4051	int tentativeStatus = problemStatus_;
4052	double changeCost;
4053	bool unflagVariables = true;
4054	bool weightsSaved = false;
4055	bool weightsSaved2 = numberIterations_ && !numberPrimalInfeasibilities_;
4056	int dontFactorizePivots = dontFactorizePivots_;
4057	if (type == 3) {
4058	type = 1;
4059	dontFactorizePivots = 1;
4060	}
4061	if (alphaAccuracy_ < 0.0 \|\| !numberPivots \|\| alphaAccuracy_ > 1.0e4 \|\| numberPivots > 20) {
4062	if (problemStatus_ > -3 \|\| numberPivots > dontFactorizePivots) {
4063	// factorize
4064	// later on we will need to recover from singularities
4065	// also we could skip if first time
4066	// save dual weights
4067	dualRowPivot_->saveWeights(this, 1);
4068	weightsSaved = true;
4069	if (type) {
4070	// is factorization okay?
4071	if (internalFactorize(1)) {
4072	// no - restore previous basis
4073	unflagVariables = false;
4074	assert (type == 1);
4075	changeMade_++; // say something changed
4076	// Keep any flagged variables
4077	int i;
4078	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4079	if (flagged(i))
4080	saveStatus_[i] \|= 64; //say flagged
4081	}
4082	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4083	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4084	numberRows_, rowActivityWork_);
4085	CoinMemcpyN(savedSolution_ ,
4086	numberColumns_, columnActivityWork_);
4087	// restore extra stuff
4088	int dummy;
4089	matrix_->generalExpanded(this, 6, dummy);
4090	// get correct bounds on all variables
4091	resetFakeBounds(1);
4092	// need to reject something
4093	char x = isColumn(sequenceOut_) ? 'C' : 'R';
4094	handler_->message(CLP_SIMPLEX_FLAG, messages_)
4095	<< x << sequenceWithin(sequenceOut_)
4096	<< CoinMessageEol;
4097	#ifdef COIN_DEVELOP
4098	printf("flag d\n");
4099	#endif
4100	setFlagged(sequenceOut_);
4101	progress_.clearBadTimes();
4102
4103	// Go to safe
4104	factorization_->pivotTolerance(0.99);
4105	forceFactorization_ = 1; // a bit drastic but ..
4106	type = 2;
4107	//assert (internalFactorize(1)==0);
4108	if (internalFactorize(1)) {
4109	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4110	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4111	numberRows_, rowActivityWork_);
4112	CoinMemcpyN(savedSolution_ ,
4113	numberColumns_, columnActivityWork_);
4114	// restore extra stuff
4115	int dummy;
4116	matrix_->generalExpanded(this, 6, dummy);
4117	// debug
4118	int returnCode = internalFactorize(1);
4119	while (returnCode) {
4120	// ouch
4121	// switch off dense
4122	int saveDense = factorization_->denseThreshold();
4123	factorization_->setDenseThreshold(0);
4124	// Go to safe
4125	factorization_->pivotTolerance(0.99);
4126	// make sure will do safe factorization
4127	pivotVariable_[0] = -1;
4128	returnCode = internalFactorize(2);
4129	factorization_->setDenseThreshold(saveDense);
4130	}
4131	// get correct bounds on all variables
4132	resetFakeBounds(1);
4133	}
4134	}
4135	}
4136	if (problemStatus_ != -4 \|\| numberPivots > 10)
4137	problemStatus_ = -3;
4138	}
4139	} else {
4140	//printf("testing with accuracy of %g and status of %d\n",alphaAccuracy_,problemStatus_);
4141	//count_alpha++;
4142	//if ((count_alpha%5000)==0)
4143	//printf("count alpha %d\n",count_alpha);
4144	}
4145	// at this stage status is -3 or -4 if looks infeasible
4146	// get primal and dual solutions
4147	#if 0
4148	{
4149	int numberTotal = numberRows_ + numberColumns_;
4150	double * saveSol = CoinCopyOfArray(solution_, numberTotal);
4151	double * saveDj = CoinCopyOfArray(dj_, numberTotal);
4152	double tolerance = type ? 1.0e-4 : 1.0e-8;
4153	// always if values pass
4154	double saveObj = objectiveValue_;
4155	double sumPrimal = sumPrimalInfeasibilities_;
4156	int numberPrimal = numberPrimalInfeasibilities_;
4157	double sumDual = sumDualInfeasibilities_;
4158	int numberDual = numberDualInfeasibilities_;
4159	gutsOfSolution(givenDuals, NULL);
4160	int j;
4161	double largestPrimal = tolerance;
4162	int iPrimal = -1;
4163	for (j = 0; j < numberTotal; j++) {
4164	double difference = solution_[j] - saveSol[j];
4165	if (fabs(difference) > largestPrimal) {
4166	iPrimal = j;
4167	largestPrimal = fabs(difference);
4168	}
4169	}
4170	double largestDual = tolerance;
4171	int iDual = -1;
4172	for (j = 0; j < numberTotal; j++) {
4173	double difference = dj_[j] - saveDj[j];
4174	if (fabs(difference) > largestDual && upper_[j] > lower_[j]) {
4175	iDual = j;
4176	largestDual = fabs(difference);
4177	}
4178	}
4179	if (!type) {
4180	if (fabs(saveObj - objectiveValue_) > 1.0e-5 \|\|
4181	numberPrimal != numberPrimalInfeasibilities_ \|\| numberPrimal != 1 \|\|
4182	fabs(sumPrimal - sumPrimalInfeasibilities_) > 1.0e-5 \|\| iPrimal >= 0 \|\|
4183	numberDual != numberDualInfeasibilities_ \|\| numberDual != 0 \|\|
4184	fabs(sumDual - sumDualInfeasibilities_) > 1.0e-5 \|\| iDual >= 0)
4185	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",
4186	type, numberIterations_, numberPivots,
4187	numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4188	largestPrimal, iPrimal,
4189	numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4190	largestDual, iDual,
4191	saveObj, objectiveValue_);
4192	} else {
4193	if (fabs(saveObj - objectiveValue_) > 1.0e-5 \|\|
4194	numberPrimalInfeasibilities_ \|\| iPrimal >= 0 \|\|
4195	numberDualInfeasibilities_ \|\| iDual >= 0)
4196	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",
4197	type, numberIterations_, numberPivots,
4198	numberPrimal, numberPrimalInfeasibilities_, sumPrimal, sumPrimalInfeasibilities_,
4199	largestPrimal, iPrimal,
4200	numberDual, numberDualInfeasibilities_, sumDual, sumDualInfeasibilities_,
4201	largestDual, iDual,
4202	saveObj, objectiveValue_);
4203	}
4204	delete [] saveSol;
4205	delete [] saveDj;
4206	}
4207	#else
4208	if (type \|\| ifValuesPass)
4209	gutsOfSolution(givenDuals, NULL);
4210	#endif
4211	// If bad accuracy treat as singular
4212	if ((largestPrimalError_ > 1.0e15 \|\| largestDualError_ > 1.0e15) && numberIterations_) {
4213	// restore previous basis
4214	unflagVariables = false;
4215	changeMade_++; // say something changed
4216	// Keep any flagged variables
4217	int i;
4218	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4219	if (flagged(i))
4220	saveStatus_[i] \|= 64; //say flagged
4221	}
4222	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4223	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4224	numberRows_, rowActivityWork_);
4225	CoinMemcpyN(savedSolution_ ,
4226	numberColumns_, columnActivityWork_);
4227	// restore extra stuff
4228	int dummy;
4229	matrix_->generalExpanded(this, 6, dummy);
4230	// get correct bounds on all variables
4231	resetFakeBounds(1);
4232	// need to reject something
4233	char x = isColumn(sequenceOut_) ? 'C' : 'R';
4234	handler_->message(CLP_SIMPLEX_FLAG, messages_)
4235	<< x << sequenceWithin(sequenceOut_)
4236	<< CoinMessageEol;
4237	#ifdef COIN_DEVELOP
4238	printf("flag e\n");
4239	#endif
4240	setFlagged(sequenceOut_);
4241	progress_.clearBadTimes();
4242
4243	// Go to safer
4244	double newTolerance = CoinMin(1.1 * factorization_->pivotTolerance(), 0.99);
4245	factorization_->pivotTolerance(newTolerance);
4246	forceFactorization_ = 1; // a bit drastic but ..
4247	if (alphaAccuracy_ != -1.0)
4248	alphaAccuracy_ = -2.0;
4249	type = 2;
4250	//assert (internalFactorize(1)==0);
4251	if (internalFactorize(1)) {
4252	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4253	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4254	numberRows_, rowActivityWork_);
4255	CoinMemcpyN(savedSolution_ ,
4256	numberColumns_, columnActivityWork_);
4257	// restore extra stuff
4258	int dummy;
4259	matrix_->generalExpanded(this, 6, dummy);
4260	// debug
4261	int returnCode = internalFactorize(1);
4262	while (returnCode) {
4263	// ouch
4264	// switch off dense
4265	int saveDense = factorization_->denseThreshold();
4266	factorization_->setDenseThreshold(0);
4267	// Go to safe
4268	factorization_->pivotTolerance(0.99);
4269	// make sure will do safe factorization
4270	pivotVariable_[0] = -1;
4271	returnCode = internalFactorize(2);
4272	factorization_->setDenseThreshold(saveDense);
4273	}
4274	// get correct bounds on all variables
4275	resetFakeBounds(1);
4276	}
4277	// get primal and dual solutions
4278	gutsOfSolution(givenDuals, NULL);
4279	} else if (goodAccuracy()) {
4280	// Can reduce tolerance
4281	double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
4282	factorization_->pivotTolerance(newTolerance);
4283	}
4284	bestObjectiveValue_ = CoinMax(bestObjectiveValue_,
4285	objectiveValue_ - bestPossibleImprovement_);
4286	bool reallyBadProblems = false;
4287	// Double check infeasibility if no action
4288	if (progress_.lastIterationNumber(0) == numberIterations_) {
4289	if (dualRowPivot_->looksOptimal()) {
4290	numberPrimalInfeasibilities_ = 0;
4291	sumPrimalInfeasibilities_ = 0.0;
4292	}
4293	#if 1
4294	} else {
4295	double thisObj = objectiveValue_ - bestPossibleImprovement_;
4296	#ifdef CLP_INVESTIGATE
4297	assert (bestPossibleImprovement_ > -1000.0 && objectiveValue_ > -1.0e100);
4298	if (bestPossibleImprovement_)
4299	printf("obj %g add in %g -> %g\n", objectiveValue_, bestPossibleImprovement_,
4300	thisObj);
4301	#endif
4302	double lastObj = progress_.lastObjective(0);
4303	#ifndef NDEBUG
4304	#ifdef COIN_DEVELOP
4305	resetFakeBounds(-1);
4306	#endif
4307	#endif
4308	#ifdef CLP_REPORT_PROGRESS
4309	ixxxxxx++;
4310	if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4311	char temp[20];
4312	sprintf(temp, "sol%d.out", ixxxxxx);
4313	printf("sol%d.out\n", ixxxxxx);
4314	FILE * fp = fopen(temp, "w");
4315	int nTotal = numberRows_ + numberColumns_;
4316	for (int i = 0; i < nTotal; i++)
4317	fprintf(fp, "%d %d %g %g %g %g %g\n",
4318	i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4319	fclose(fp);
4320	}
4321	#endif
4322	if(!ifValuesPass && firstFree_ < 0) {
4323	double testTol = 5.0e-3;
4324	if (progress_.timesFlagged() > 10) {
4325	testTol *= pow(2.0, progress_.timesFlagged() - 8);
4326	} else if (progress_.timesFlagged() > 5) {
4327	testTol *= 5.0;
4328	}
4329	if (lastObj > thisObj +
4330	testTol*(fabs(thisObj) + fabs(lastObj)) + testTol) {
4331	int maxFactor = factorization_->maximumPivots();
4332	if ((specialOptions_ & 1048576) == 0) {
4333	if (progress_.timesFlagged() > 10)
4334	progress_.incrementReallyBadTimes();
4335	if (maxFactor > 10 - 9) {
4336	#ifdef COIN_DEVELOP
4337	printf("lastobj %g thisobj %g\n", lastObj, thisObj);
4338	#endif
4339	//if (forceFactorization_<0)
4340	//forceFactorization_= maxFactor;
4341	//forceFactorization_ = CoinMax(1,(forceFactorization_>>1));
4342	if ((progressFlag_ & 4) == 0 && lastObj < thisObj + 1.0e4 &&
4343	largestPrimalError_ < 1.0e2) {
4344	// Just save costs
4345	// save extra copy of cost_
4346	int nTotal = numberRows_ + numberColumns_;
4347	double * temp = new double [2*nTotal];
4348	memcpy(temp, cost_, nTotal * sizeof(double));
4349	memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4350	delete [] cost_;
4351	cost_ = temp;
4352	objectiveWork_ = cost_;
4353	rowObjectiveWork_ = cost_ + numberColumns_;
4354	progressFlag_ \|= 4;
4355	} else {
4356	forceFactorization_ = 1;
4357	#ifdef COIN_DEVELOP
4358	printf("Reducing factorization frequency - bad backwards\n");
4359	#endif
4360	#if 1
4361	unflagVariables = false;
4362	changeMade_++; // say something changed
4363	int nTotal = numberRows_ + numberColumns_;
4364	CoinMemcpyN(saveStatus_, nTotal, status_);
4365	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4366	numberRows_, rowActivityWork_);
4367	CoinMemcpyN(savedSolution_ ,
4368	numberColumns_, columnActivityWork_);
4369	if ((progressFlag_ & 4) == 0) {
4370	// save extra copy of cost_
4371	double * temp = new double [2*nTotal];
4372	memcpy(temp, cost_, nTotal * sizeof(double));
4373	memcpy(temp + nTotal, cost_, nTotal * sizeof(double));
4374	delete [] cost_;
4375	cost_ = temp;
4376	objectiveWork_ = cost_;
4377	rowObjectiveWork_ = cost_ + numberColumns_;
4378	progressFlag_ \|= 4;
4379	} else {
4380	memcpy(cost_, cost_ + nTotal, nTotal * sizeof(double));
4381	}
4382	// restore extra stuff
4383	int dummy;
4384	matrix_->generalExpanded(this, 6, dummy);
4385	double pivotTolerance = factorization_->pivotTolerance();
4386	if(pivotTolerance < 0.2)
4387	factorization_->pivotTolerance(0.2);
4388	else if(progress_.timesFlagged() > 2)
4389	factorization_->pivotTolerance(CoinMin(pivotTolerance * 1.1, 0.99));
4390	if (alphaAccuracy_ != -1.0)
4391	alphaAccuracy_ = -2.0;
4392	if (internalFactorize(1)) {
4393	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
4394	CoinMemcpyN(savedSolution_ + numberColumns_ ,
4395	numberRows_, rowActivityWork_);
4396	CoinMemcpyN(savedSolution_ ,
4397	numberColumns_, columnActivityWork_);
4398	// restore extra stuff
4399	int dummy;
4400	matrix_->generalExpanded(this, 6, dummy);
4401	// debug
4402	int returnCode = internalFactorize(1);
4403	while (returnCode) {
4404	// ouch
4405	// switch off dense
4406	int saveDense = factorization_->denseThreshold();
4407	factorization_->setDenseThreshold(0);
4408	// Go to safe
4409	factorization_->pivotTolerance(0.99);
4410	// make sure will do safe factorization
4411	pivotVariable_[0] = -1;
4412	returnCode = internalFactorize(2);
4413	factorization_->setDenseThreshold(saveDense);
4414	}
4415	}
4416	resetFakeBounds(0);
4417	type = 2; // so will restore weights
4418	// get primal and dual solutions
4419	gutsOfSolution(givenDuals, NULL);
4420	if (numberPivots < 2) {
4421	// need to reject something
4422	char x = isColumn(sequenceOut_) ? 'C' : 'R';
4423	handler_->message(CLP_SIMPLEX_FLAG, messages_)
4424	<< x << sequenceWithin(sequenceOut_)
4425	<< CoinMessageEol;
4426	#ifdef COIN_DEVELOP
4427	printf("flag d\n");
4428	#endif
4429	setFlagged(sequenceOut_);
4430	progress_.clearBadTimes();
4431	progress_.incrementTimesFlagged();
4432	}
4433	if (numberPivots < 10)
4434	reallyBadProblems = true;
4435	#ifdef COIN_DEVELOP
4436	printf("obj now %g\n", objectiveValue_);
4437	#endif
4438	progress_.modifyObjective(objectiveValue_
4439	- bestPossibleImprovement_);
4440	#endif
4441	}
4442	}
4443	} else {
4444	// in fast dual give up
4445	#ifdef COIN_DEVELOP
4446	printf("In fast dual?\n");
4447	#endif
4448	problemStatus_ = 3;
4449	}
4450	} else if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-3) {
4451	numberTimesOptimal_ = 0;
4452	}
4453	}
4454	#endif
4455	}
4456	// Up tolerance if looks a bit odd
4457	if (numberIterations_ > CoinMax(1000, numberRows_ >> 4) && (specialOptions_ & 64) != 0) {
4458	if (sumPrimalInfeasibilities_ && sumPrimalInfeasibilities_ < 1.0e5) {
4459	int backIteration = progress_.lastIterationNumber(CLP_PROGRESS - 1);
4460	if (backIteration > 0 && numberIterations_ - backIteration < 9 * CLP_PROGRESS) {
4461	if (factorization_->pivotTolerance() < 0.9) {
4462	// up tolerance
4463	factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
4464	//printf("tol now %g\n",factorization_->pivotTolerance());
4465	progress_.clearIterationNumbers();
4466	}
4467	}
4468	}
4469	}
4470	// Check if looping
4471	int loop;
4472	if (!givenDuals && type != 2)
4473	loop = progress_.looping();
4474	else
4475	loop = -1;
4476	if (progress_.reallyBadTimes() > 10) {
4477	problemStatus_ = 10; // instead - try other algorithm
4478	#if COIN_DEVELOP>2
4479	printf("returning at %d\n", __LINE__);
4480	#endif
4481	}
4482	int situationChanged = 0;
4483	if (loop >= 0) {
4484	problemStatus_ = loop; //exit if in loop
4485	if (!problemStatus_) {
4486	// declaring victory
4487	numberPrimalInfeasibilities_ = 0;
4488	sumPrimalInfeasibilities_ = 0.0;
4489	} else {
4490	problemStatus_ = 10; // instead - try other algorithm
4491	#if COIN_DEVELOP>2
4492	printf("returning at %d\n", __LINE__);
4493	#endif
4494	}
4495	return;
4496	} else if (loop < -1) {
4497	// something may have changed
4498	gutsOfSolution(NULL, NULL);
4499	situationChanged = 1;
4500	}
4501	// really for free variables in
4502	if((progressFlag_ & 2) != 0) {
4503	situationChanged = 2;
4504	}
4505	progressFlag_ &= (~3); //reset progress flag
4506	if ((progressFlag_ & 4) != 0) {
4507	// save copy of cost_
4508	int nTotal = numberRows_ + numberColumns_;
4509	memcpy(cost_ + nTotal, cost_, nTotal * sizeof(double));
4510	}
4511	/*if (!numberIterations_&&sumDualInfeasibilities_)
4512	printf("OBJ %g sumPinf %g sumDinf %g\n",
4513	objectiveValue(),sumPrimalInfeasibilities_,
4514	sumDualInfeasibilities_);*/
4515	// mark as having gone optimal if looks like it
4516	if (!numberPrimalInfeasibilities_&&
4517	!numberDualInfeasibilities_)
4518	progressFlag_ \|= 8;
4519	if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
4520	handler_->message(CLP_SIMPLEX_STATUS, messages_)
4521	<< numberIterations_ << objectiveValue();
4522	handler_->printing(sumPrimalInfeasibilities_ > 0.0)
4523	<< sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
4524	handler_->printing(sumDualInfeasibilities_ > 0.0)
4525	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
4526	handler_->printing(numberDualInfeasibilitiesWithoutFree_
4527	< numberDualInfeasibilities_)
4528	<< numberDualInfeasibilitiesWithoutFree_;
4529	handler_->message() << CoinMessageEol;
4530	}
4531	#if 0
4532	printf("IT %d %g %g(%d) %g(%d)\n",
4533	numberIterations_, objectiveValue(),
4534	sumPrimalInfeasibilities_, numberPrimalInfeasibilities_,
4535	sumDualInfeasibilities_, numberDualInfeasibilities_);
4536	#endif
4537	double approximateObjective = objectiveValue_;
4538	#ifdef CLP_REPORT_PROGRESS
4539	if (ixxxxxx >= ixxyyyy - 4 && ixxxxxx <= ixxyyyy) {
4540	char temp[20];
4541	sprintf(temp, "x_sol%d.out", ixxxxxx);
4542	FILE * fp = fopen(temp, "w");
4543	int nTotal = numberRows_ + numberColumns_;
4544	for (int i = 0; i < nTotal; i++)
4545	fprintf(fp, "%d %d %g %g %g %g %g\n",
4546	i, status_[i], lower_[i], solution_[i], upper_[i], cost_[i], dj_[i]);
4547	fclose(fp);
4548	if (ixxxxxx == ixxyyyy)
4549	exit(6);
4550	}
4551	#endif
4552	realDualInfeasibilities = sumDualInfeasibilities_;
4553	double saveTolerance = dualTolerance_;
4554	// If we need to carry on cleaning variables
4555	if (!numberPrimalInfeasibilities_ && (specialOptions_ & 1024) != 0 && CLEAN_FIXED) {
4556	for (int iRow = 0; iRow < numberRows_; iRow++) {
4557	int iPivot = pivotVariable_[iRow];
4558	if (!flagged(iPivot) && pivoted(iPivot)) {
4559	// carry on
4560	numberPrimalInfeasibilities_ = -1;
4561	sumOfRelaxedPrimalInfeasibilities_ = 1.0;
4562	sumPrimalInfeasibilities_ = 1.0;
4563	break;
4564	}
4565	}
4566	}
4567	/* If we are primal feasible and any dual infeasibilities are on
4568	free variables then it is better to go to primal */
4569	if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
4570	numberDualInfeasibilities_)
4571	problemStatus_ = 10;
4572	// dual bound coming in
4573	//double saveDualBound = dualBound_;
4574	bool needCleanFake = false;
4575	while (problemStatus_ <= -3) {
4576	int cleanDuals = 0;
4577	if (situationChanged != 0)
4578	cleanDuals = 1;
4579	int numberChangedBounds = 0;
4580	int doOriginalTolerance = 0;
4581	if ( lastCleaned == numberIterations_)
4582	doOriginalTolerance = 1;
4583	// check optimal
4584	// give code benefit of doubt
4585	if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
4586	sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
4587	// say optimal (with these bounds etc)
4588	numberDualInfeasibilities_ = 0;
4589	sumDualInfeasibilities_ = 0.0;
4590	numberPrimalInfeasibilities_ = 0;
4591	sumPrimalInfeasibilities_ = 0.0;
4592	}
4593	//if (dualFeasible()\|\|problemStatus_==-4\|\|(primalFeasible()&&!numberDualInfeasibilitiesWithoutFree_)) {
4594	if (dualFeasible() \|\| problemStatus_ == -4) {
4595	progress_.modifyObjective(objectiveValue_
4596	- bestPossibleImprovement_);
4597	#ifdef COIN_DEVELOP
4598	if (sumDualInfeasibilities_ \|\| bestPossibleImprovement_)
4599	printf("improve %g dualinf %g -> %g\n",
4600	bestPossibleImprovement_, sumDualInfeasibilities_,
4601	sumDualInfeasibilities_ * dualBound_);
4602	#endif
4603	// see if cutoff reached
4604	double limit = 0.0;
4605	getDblParam(ClpDualObjectiveLimit, limit);
4606	#if 0
4607	if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
4608	limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberAtFakeBound()) {
4609	//looks infeasible on objective
4610	if (perturbation_ == 101) {
4611	cleanDuals = 1;
4612	// Save costs
4613	int numberTotal = numberRows_ + numberColumns_;
4614	double * saveCost = CoinCopyOfArray(cost_, numberTotal);
4615	// make sure fake bounds are back
4616	changeBounds(1, NULL, changeCost);
4617	createRim4(false);
4618	// make sure duals are current
4619	computeDuals(givenDuals);
4620	checkDualSolution();
4621	if(objectiveValue()*optimizationDirection_ >
4622	limit + 1.0e-7 + 1.0e-8 * fabs(limit) && !numberDualInfeasibilities_) {
4623	perturbation_ = 102; // stop any perturbations
4624	printf("cutoff test succeeded\n");
4625	} else {
4626	printf("cutoff test failed\n");
4627	// put back
4628	memcpy(cost_, saveCost, numberTotal * sizeof(double));
4629	// make sure duals are current
4630	computeDuals(givenDuals);
4631	checkDualSolution();
4632	progress_.modifyObjective(-COIN_DBL_MAX);
4633	problemStatus_ = -1;
4634	}
4635	delete [] saveCost;
4636	}
4637	}
4638	#endif
4639	if (primalFeasible() && !givenDuals) {
4640	// may be optimal - or may be bounds are wrong
4641	handler_->message(CLP_DUAL_BOUNDS, messages_)
4642	<< dualBound_
4643	<< CoinMessageEol;
4644	// save solution in case unbounded
4645	double * saveColumnSolution = NULL;
4646	double * saveRowSolution = NULL;
4647	bool inCbc = (specialOptions_ & (0x01000000 \| 16384)) != 0;
4648	if (!inCbc) {
4649	saveColumnSolution = CoinCopyOfArray(columnActivityWork_, numberColumns_);
4650	saveRowSolution = CoinCopyOfArray(rowActivityWork_, numberRows_);
4651	}
4652	numberChangedBounds = changeBounds(0, rowArray_[3], changeCost);
4653	if (numberChangedBounds <= 0 && !numberDualInfeasibilities_) {
4654	//looks optimal - do we need to reset tolerance
4655	if (perturbation_ == 101) {
4656	perturbation_ = 102; // stop any perturbations
4657	cleanDuals = 1;
4658	// make sure fake bounds are back
4659	//computeObjectiveValue();
4660	changeBounds(1, NULL, changeCost);
4661	//computeObjectiveValue();
4662	createRim4(false);
4663	// make sure duals are current
4664	computeDuals(givenDuals);
4665	checkDualSolution();
4666	progress_.modifyObjective(-COIN_DBL_MAX);
4667	#define DUAL_TRY_FASTER
4668	#ifdef DUAL_TRY_FASTER
4669	if (numberDualInfeasibilities_) {
4670	#endif
4671	numberChanged_ = 1; // force something to happen
4672	lastCleaned = numberIterations_ - 1;
4673	#ifdef DUAL_TRY_FASTER
4674	} else {
4675	//double value = objectiveValue_;
4676	computeObjectiveValue(true);
4677	//printf("old %g new %g\n",value,objectiveValue_);
4678	//numberChanged_=1;
4679	}
4680	#endif
4681	}
4682	if (lastCleaned < numberIterations_ && numberTimesOptimal_ < 4 &&
4683	(numberChanged_ \|\| (specialOptions_ & 4096) == 0)) {
4684	#if CLP_CAN_HAVE_ZERO_OBJ
4685	if ((specialOptions_&2097152)==0) {
4686	#endif
4687	doOriginalTolerance = 2;
4688	numberTimesOptimal_++;
4689	changeMade_++; // say something changed
4690	if (numberTimesOptimal_ == 1) {
4691	dualTolerance_ = dblParam_[ClpDualTolerance];
4692	} else {
4693	if (numberTimesOptimal_ == 2) {
4694	// better to have small tolerance even if slower
4695	factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
4696	}
4697	dualTolerance_ = dblParam_[ClpDualTolerance];
4698	dualTolerance_ *= pow(2.0, numberTimesOptimal_ - 1);
4699	}
4700	cleanDuals = 2; // If nothing changed optimal else primal
4701	#if CLP_CAN_HAVE_ZERO_OBJ
4702	} else {
4703	// no cost - skip checks
4704	problemStatus_=0;
4705	}
4706	#endif
4707	} else {
4708	problemStatus_ = 0; // optimal
4709	if (lastCleaned < numberIterations_ && numberChanged_) {
4710	handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
4711	<< CoinMessageEol;
4712	}
4713	}
4714	} else {
4715	cleanDuals = 1;
4716	if (doOriginalTolerance == 1) {
4717	// check unbounded
4718	// find a variable with bad dj
4719	int iSequence;
4720	int iChosen = -1;
4721	if (!inCbc) {
4722	double largest = 100.0 * primalTolerance_;
4723	for (iSequence = 0; iSequence < numberRows_ + numberColumns_;
4724	iSequence++) {
4725	double djValue = dj_[iSequence];
4726	double originalLo = originalLower(iSequence);
4727	double originalUp = originalUpper(iSequence);
4728	if (fabs(djValue) > fabs(largest)) {
4729	if (getStatus(iSequence) != basic) {
4730	if (djValue > 0 && originalLo < -1.0e20) {
4731	if (djValue > fabs(largest)) {
4732	largest = djValue;
4733	iChosen = iSequence;
4734	}
4735	} else if (djValue < 0 && originalUp > 1.0e20) {
4736	if (-djValue > fabs(largest)) {
4737	largest = djValue;
4738	iChosen = iSequence;
4739	}
4740	}
4741	}
4742	}
4743	}
4744	}
4745	if (iChosen >= 0) {
4746	int iSave = sequenceIn_;
4747	sequenceIn_ = iChosen;
4748	unpack(rowArray_[1]);
4749	sequenceIn_ = iSave;
4750	// if dual infeasibilities then must be free vector so add in dual
4751	if (numberDualInfeasibilities_) {
4752	if (fabs(changeCost) > 1.0e-5)
4753	COIN_DETAIL_PRINT(printf("Odd free/unbounded combo\n"));
4754	changeCost += cost_[iChosen];
4755	}
4756	problemStatus_ = checkUnbounded(rowArray_[1], rowArray_[0],
4757	changeCost);
4758	rowArray_[1]->clear();
4759	} else {
4760	problemStatus_ = -3;
4761	}
4762	if (problemStatus_ == 2 && perturbation_ == 101) {
4763	perturbation_ = 102; // stop any perturbations
4764	cleanDuals = 1;
4765	createRim4(false);
4766	progress_.modifyObjective(-COIN_DBL_MAX);
4767	problemStatus_ = -1;
4768	}
4769	if (problemStatus_ == 2) {
4770	// it is unbounded - restore solution
4771	// but first add in changes to non-basic
4772	int iColumn;
4773	double * original = columnArray_[0]->denseVector();
4774	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4775	if(getColumnStatus(iColumn) != basic)
4776	ray_[iColumn] +=
4777	saveColumnSolution[iColumn] - original[iColumn];
4778	columnActivityWork_[iColumn] = original[iColumn];
4779	}
4780	CoinMemcpyN(saveRowSolution, numberRows_,
4781	rowActivityWork_);
4782	}
4783	} else {
4784	doOriginalTolerance = 2;
4785	rowArray_[0]->clear();
4786	}
4787	}
4788	delete [] saveColumnSolution;
4789	delete [] saveRowSolution;
4790	}
4791	if (problemStatus_ == -4 \|\| problemStatus_ == -5) {
4792	// may be infeasible - or may be bounds are wrong
4793	numberChangedBounds = changeBounds(0, NULL, changeCost);
4794	needCleanFake = true;
4795	/* Should this be here as makes no difference to being feasible.
4796	But seems to make a difference to run times. */
4797	if (perturbation_ == 101 && 0) {
4798	perturbation_ = 102; // stop any perturbations
4799	cleanDuals = 1;
4800	numberChangedBounds = 1;
4801	// make sure fake bounds are back
4802	changeBounds(1, NULL, changeCost);
4803	needCleanFake = true;
4804	createRim4(false);
4805	progress_.modifyObjective(-COIN_DBL_MAX);
4806	}
4807	if ((numberChangedBounds <= 0 \|\| dualBound_ > 1.0e20 \|\|
4808	(largestPrimalError_ > 1.0 && dualBound_ > 1.0e17)) &&
4809	(numberPivots < 4 \|\| sumPrimalInfeasibilities_ > 1.0e-6)) {
4810	problemStatus_ = 1; // infeasible
4811	if (perturbation_ == 101) {
4812	perturbation_ = 102; // stop any perturbations
4813	//cleanDuals=1;
4814	//numberChangedBounds=1;
4815	//createRim4(false);
4816	}
4817	} else {
4818	problemStatus_ = -1; //iterate
4819	cleanDuals = 1;
4820	if (numberChangedBounds <= 0)
4821	doOriginalTolerance = 2;
4822	// and delete ray which has been created
4823	delete [] ray_;
4824	ray_ = NULL;
4825	}
4826
4827	}
4828	} else {
4829	cleanDuals = 1;
4830	}
4831	if (problemStatus_ < 0) {
4832	if (doOriginalTolerance == 2) {
4833	// put back original tolerance
4834	lastCleaned = numberIterations_;
4835	numberChanged_ = 0; // Number of variables with changed costs
4836	handler_->message(CLP_DUAL_ORIGINAL, messages_)
4837	<< CoinMessageEol;
4838	perturbation_ = 102; // stop any perturbations
4839	#if 0
4840	double * xcost = new double[numberRows_+numberColumns_];
4841	double * xlower = new double[numberRows_+numberColumns_];
4842	double * xupper = new double[numberRows_+numberColumns_];
4843	double * xdj = new double[numberRows_+numberColumns_];
4844	double * xsolution = new double[numberRows_+numberColumns_];
4845	CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4846	CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4847	CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4848	CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4849	CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4850	#endif
4851	createRim4(false);
4852	progress_.modifyObjective(-COIN_DBL_MAX);
4853	// make sure duals are current
4854	computeDuals(givenDuals);
4855	checkDualSolution();
4856	#if 0
4857	int i;
4858	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4859	if (cost_[i] != xcost[i])
4860	printf("** %d old cost %g new %g sol %g\n",
4861	i, xcost[i], cost_[i], solution_[i]);
4862	if (lower_[i] != xlower[i])
4863	printf("** %d old lower %g new %g sol %g\n",
4864	i, xlower[i], lower_[i], solution_[i]);
4865	if (upper_[i] != xupper[i])
4866	printf("** %d old upper %g new %g sol %g\n",
4867	i, xupper[i], upper_[i], solution_[i]);
4868	if (dj_[i] != xdj[i])
4869	printf("** %d old dj %g new %g sol %g\n",
4870	i, xdj[i], dj_[i], solution_[i]);
4871	if (solution_[i] != xsolution[i])
4872	printf("** %d old solution %g new %g sol %g\n",
4873	i, xsolution[i], solution_[i], solution_[i]);
4874	}
4875	//delete [] xcost;
4876	//delete [] xupper;
4877	//delete [] xlower;
4878	//delete [] xdj;
4879	//delete [] xsolution;
4880	#endif
4881	// put back bounds as they were if was optimal
4882	if (doOriginalTolerance == 2 && cleanDuals != 2) {
4883	changeMade_++; // say something changed
4884	/* We may have already changed some bounds in this function
4885	so save numberFake_ and add in.
4886
4887	Worst that can happen is that we waste a bit of time - but it must be finite.
4888	*/
4889	//int saveNumberFake = numberFake_;
4890	//resetFakeBounds(-1);
4891	changeBounds(3, NULL, changeCost);
4892	needCleanFake = true;
4893	//numberFake_ += saveNumberFake;
4894	//resetFakeBounds(-1);
4895	cleanDuals = 2;
4896	//cleanDuals=1;
4897	}
4898	#if 0
4899	//int i;
4900	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4901	if (cost_[i] != xcost[i])
4902	printf("** %d old cost %g new %g sol %g\n",
4903	i, xcost[i], cost_[i], solution_[i]);
4904	if (lower_[i] != xlower[i])
4905	printf("** %d old lower %g new %g sol %g\n",
4906	i, xlower[i], lower_[i], solution_[i]);
4907	if (upper_[i] != xupper[i])
4908	printf("** %d old upper %g new %g sol %g\n",
4909	i, xupper[i], upper_[i], solution_[i]);
4910	if (dj_[i] != xdj[i])
4911	printf("** %d old dj %g new %g sol %g\n",
4912	i, xdj[i], dj_[i], solution_[i]);
4913	if (solution_[i] != xsolution[i])
4914	printf("** %d old solution %g new %g sol %g\n",
4915	i, xsolution[i], solution_[i], solution_[i]);
4916	}
4917	delete [] xcost;
4918	delete [] xupper;
4919	delete [] xlower;
4920	delete [] xdj;
4921	delete [] xsolution;
4922	#endif
4923	}
4924	if (cleanDuals == 1 \|\| (cleanDuals == 2 && !numberDualInfeasibilities_)) {
4925	// make sure dual feasible
4926	// look at all rows and columns
4927	rowArray_[0]->clear();
4928	columnArray_[0]->clear();
4929	double objectiveChange = 0.0;
4930	double savePrimalInfeasibilities = sumPrimalInfeasibilities_;
4931	if (!numberIterations_) {
4932	int nTotal = numberRows_ + numberColumns_;
4933	if (arraysNotCreated) {
4934	// create save arrays
4935	delete [] saveStatus_;
4936	delete [] savedSolution_;
4937	saveStatus_ = new unsigned char [nTotal];
4938	savedSolution_ = new double [nTotal];
4939	arraysNotCreated = false;
4940	}
4941	// save arrays
4942	CoinMemcpyN(status_, nTotal, saveStatus_);
4943	CoinMemcpyN(rowActivityWork_,
4944	numberRows_, savedSolution_ + numberColumns_);
4945	CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
4946	}
4947	#if 0
4948	double * xcost = new double[numberRows_+numberColumns_];
4949	double * xlower = new double[numberRows_+numberColumns_];
4950	double * xupper = new double[numberRows_+numberColumns_];
4951	double * xdj = new double[numberRows_+numberColumns_];
4952	double * xsolution = new double[numberRows_+numberColumns_];
4953	CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
4954	CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
4955	CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
4956	CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
4957	CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
4958	#endif
4959	if (givenDuals)
4960	dualTolerance_ = 1.0e50;
4961	#if CLP_CAN_HAVE_ZERO_OBJ>1
4962	if ((specialOptions_&2097152)==0) {
4963	#endif
4964	updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
4965	0.0, objectiveChange, true);
4966	#if CLP_CAN_HAVE_ZERO_OBJ>1
4967	} else {
4968	rowArray_[0]->clear();
4969	rowArray_[1]->clear();
4970	columnArray_[0]->clear();
4971	}
4972	#endif
4973	dualTolerance_ = saveTolerance;
4974	#if 0
4975	int i;
4976	for (i = 0; i < numberRows_ + numberColumns_; i++) {
4977	if (cost_[i] != xcost[i])
4978	printf("** %d old cost %g new %g sol %g\n",
4979	i, xcost[i], cost_[i], solution_[i]);
4980	if (lower_[i] != xlower[i])
4981	printf("** %d old lower %g new %g sol %g\n",
4982	i, xlower[i], lower_[i], solution_[i]);
4983	if (upper_[i] != xupper[i])
4984	printf("** %d old upper %g new %g sol %g\n",
4985	i, xupper[i], upper_[i], solution_[i]);
4986	if (dj_[i] != xdj[i])
4987	printf("** %d old dj %g new %g sol %g\n",
4988	i, xdj[i], dj_[i], solution_[i]);
4989	if (solution_[i] != xsolution[i])
4990	printf("** %d old solution %g new %g sol %g\n",
4991	i, xsolution[i], solution_[i], solution_[i]);
4992	}
4993	delete [] xcost;
4994	delete [] xupper;
4995	delete [] xlower;
4996	delete [] xdj;
4997	delete [] xsolution;
4998	#endif
4999	// for now - recompute all
5000	gutsOfSolution(NULL, NULL);
5001	if (givenDuals)
5002	dualTolerance_ = 1.0e50;
5003	#if CLP_CAN_HAVE_ZERO_OBJ>1
5004	if ((specialOptions_&2097152)==0) {
5005	#endif
5006	updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
5007	0.0, objectiveChange, true);
5008	#if CLP_CAN_HAVE_ZERO_OBJ>1
5009	} else {
5010	rowArray_[0]->clear();
5011	rowArray_[1]->clear();
5012	columnArray_[0]->clear();
5013	}
5014	#endif
5015	dualTolerance_ = saveTolerance;
5016	if (!numberIterations_ && sumPrimalInfeasibilities_ >
5017	1.0e5*(savePrimalInfeasibilities+1.0e3) &&
5018	(moreSpecialOptions_ & 256) == 0) {
5019	// Use primal
5020	int nTotal = numberRows_ + numberColumns_;
5021	CoinMemcpyN(saveStatus_, nTotal, status_);
5022	CoinMemcpyN(savedSolution_ + numberColumns_ ,
5023	numberRows_, rowActivityWork_);
5024	CoinMemcpyN(savedSolution_ ,
5025	numberColumns_, columnActivityWork_);
5026	problemStatus_ = 10;
5027	situationChanged = 0;
5028	}
5029	//assert(numberDualInfeasibilitiesWithoutFree_==0);
5030	if (numberDualInfeasibilities_) {
5031	if ((numberPrimalInfeasibilities_ \|\| numberPivots)
5032	&& problemStatus_!=10) {
5033	problemStatus_ = -1; // carry on as normal
5034	} else {
5035	problemStatus_ = 10; // try primal
5036	#if COIN_DEVELOP>1
5037	printf("returning at %d\n", __LINE__);
5038	#endif
5039	}
5040	} else if (situationChanged == 2) {
5041	problemStatus_ = -1; // carry on as normal
5042	// need to reset bounds
5043	changeBounds(3, NULL, changeCost);
5044	}
5045	situationChanged = 0;
5046	} else {
5047	// iterate
5048	if (cleanDuals != 2) {
5049	problemStatus_ = -1;
5050	} else {
5051	problemStatus_ = 10; // try primal
5052	#if COIN_DEVELOP>2
5053	printf("returning at %d\n", __LINE__);
5054	#endif
5055	}
5056	}
5057	}
5058	}
5059	// unflag all variables (we may want to wait a bit?)
5060	if ((tentativeStatus != -2 && tentativeStatus != -1) && unflagVariables) {
5061	int iRow;
5062	int numberFlagged = 0;
5063	for (iRow = 0; iRow < numberRows_; iRow++) {
5064	int iPivot = pivotVariable_[iRow];
5065	if (flagged(iPivot)) {
5066	numberFlagged++;
5067	clearFlagged(iPivot);
5068	}
5069	}
5070	#ifdef COIN_DEVELOP
5071	if (numberFlagged) {
5072	printf("unflagging %d variables - tentativeStatus %d probStat %d ninf %d nopt %d\n", numberFlagged, tentativeStatus,
5073	problemStatus_, numberPrimalInfeasibilities_,
5074	numberTimesOptimal_);
5075	}
5076	#endif
5077	unflagVariables = numberFlagged > 0;
5078	if (numberFlagged && !numberPivots) {
5079	/* looks like trouble as we have not done any iterations.
5080	Try changing pivot tolerance then give it a few goes and give up */
5081	if (factorization_->pivotTolerance() < 0.9) {
5082	factorization_->pivotTolerance(0.99);
5083	problemStatus_ = -1;
5084	} else if (numberTimesOptimal_ < 3) {
5085	numberTimesOptimal_++;
5086	problemStatus_ = -1;
5087	} else {
5088	unflagVariables = false;
5089	//secondaryStatus_ = 1; // and say probably infeasible
5090	if ((moreSpecialOptions_ & 256) == 0) {
5091	// try primal
5092	problemStatus_ = 10;
5093	} else {
5094	// almost certainly infeasible
5095	problemStatus_ = 1;
5096	}
5097	#if COIN_DEVELOP>1
5098	printf("returning at %d\n", __LINE__);
5099	#endif
5100	}
5101	}
5102	}
5103	if (problemStatus_ < 0) {
5104	if (needCleanFake) {
5105	double dummyChangeCost = 0.0;
5106	changeBounds(3, NULL, dummyChangeCost);
5107	}
5108	#if 0
5109	if (objectiveValue_ < lastObjectiveValue_ - 1.0e-8 *
5110	CoinMax(fabs(objectivevalue_), fabs(lastObjectiveValue_))) {
5111	} else {
5112	lastObjectiveValue_ = objectiveValue_;
5113	}
5114	#endif
5115	if (type == 0 \|\| type == 1) {
5116	if (!type && arraysNotCreated) {
5117	// create save arrays
5118	delete [] saveStatus_;
5119	delete [] savedSolution_;
5120	saveStatus_ = new unsigned char [numberRows_+numberColumns_];
5121	savedSolution_ = new double [numberRows_+numberColumns_];
5122	}
5123	// save arrays
5124	CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
5125	CoinMemcpyN(rowActivityWork_,
5126	numberRows_, savedSolution_ + numberColumns_);
5127	CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
5128	// save extra stuff
5129	int dummy;
5130	matrix_->generalExpanded(this, 5, dummy);
5131	}
5132	if (weightsSaved) {
5133	// restore weights (if saved) - also recompute infeasibility list
5134	if (!reallyBadProblems && (largestPrimalError_ < 100.0 \|\| numberPivots > 10)) {
5135	if (tentativeStatus > -3)
5136	dualRowPivot_->saveWeights(this, (type < 2) ? 2 : 4);
5137	else
5138	dualRowPivot_->saveWeights(this, 3);
5139	} else {
5140	// reset weights or scale back
5141	dualRowPivot_->saveWeights(this, 6);
5142	}
5143	} else if (weightsSaved2 && numberPrimalInfeasibilities_) {
5144	dualRowPivot_->saveWeights(this, 3);
5145	}
5146	}
5147	// see if cutoff reached
5148	double limit = 0.0;
5149	getDblParam(ClpDualObjectiveLimit, limit);
5150	#if 0
5151	if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5152	limit + 100.0) {
5153	printf("lim %g obj %g %g - wo perturb %g sum dual %g\n",
5154	limit, objectiveValue_, objectiveValue(), computeInternalObjectiveValue(), sumDualInfeasibilities_);
5155	}
5156	#endif
5157	if(fabs(limit) < 1.0e30 && objectiveValue()*optimizationDirection_ >
5158	limit && !numberAtFakeBound()) {
5159	bool looksInfeasible = !numberDualInfeasibilities_;
5160	if (objectiveValue()optimizationDirection_ > limit + fabs(0.1 limit) + 1.0e2 * sumDualInfeasibilities_ + 1.0e4 &&
5161	sumDualInfeasibilities_ < largestDualError_ && numberIterations_ > 0.5 * numberRows_ + 1000)
5162	looksInfeasible = true;
5163	if (looksInfeasible) {
5164	// Even if not perturbed internal costs may have changed
5165	// be careful
5166	if (true \|\| numberIterations_) {
5167	if(computeInternalObjectiveValue() > limit) {
5168	problemStatus_ = 1;
5169	secondaryStatus_ = 1; // and say was on cutoff
5170	}
5171	} else {
5172	problemStatus_ = 1;
5173	secondaryStatus_ = 1; // and say was on cutoff
5174	}
5175	}
5176	}
5177	// If we are in trouble and in branch and bound give up
5178	if ((specialOptions_ & 1024) != 0) {
5179	int looksBad = 0;
5180	if (largestPrimalError_ * largestDualError_ > 1.0e2) {
5181	looksBad = 1;
5182	} else if (largestPrimalError_ > 1.0e-2
5183	&& objectiveValue_ > CoinMin(1.0e15, 1.0e3 * limit)) {
5184	looksBad = 2;
5185	}
5186	if (looksBad) {
5187	if (factorization_->pivotTolerance() < 0.9) {
5188	// up tolerance
5189	factorization_->pivotTolerance(CoinMin(factorization_->pivotTolerance() * 1.05 + 0.02, 0.91));
5190	} else if (numberIterations_ > 10000) {
5191	if (handler_->logLevel() > 2)
5192	printf("bad dual - saying infeasible %d\n", looksBad);
5193	problemStatus_ = 1;
5194	secondaryStatus_ = 1; // and say was on cutoff
5195	} else if (largestPrimalError_ > 1.0e5) {
5196	{
5197	int iBigB = -1;
5198	double bigB = 0.0;
5199	int iBigN = -1;
5200	double bigN = 0.0;
5201	for (int i = 0; i < numberRows_ + numberColumns_; i++) {
5202	double value = fabs(solution_[i]);
5203	if (getStatus(i) == basic) {
5204	if (value > bigB) {
5205	bigB = value;
5206	iBigB = i;
5207	}
5208	} else {
5209	if (value > bigN) {
5210	bigN = value;
5211	iBigN = i;
5212	}
5213	}
5214	}
5215	#ifdef CLP_INVESTIGATE
5216	if (bigB > 1.0e8 \|\| bigN > 1.0e8) {
5217	if (handler_->logLevel() > 0)
5218	printf("it %d - basic %d %g, nonbasic %d %g\n",
5219	numberIterations_, iBigB, bigB, iBigN, bigN);
5220	}
5221	#endif
5222	}
5223	#if COIN_DEVELOP!=2
5224	if (handler_->logLevel() > 2)
5225	#endif
5226	printf("bad dual - going to primal %d %g\n", looksBad, largestPrimalError_);
5227	allSlackBasis(true);
5228	problemStatus_ = 10;
5229	}
5230	}
5231	}
5232	if (problemStatus_ < 0 && !changeMade_) {
5233	problemStatus_ = 4; // unknown
5234	}
5235	lastGoodIteration_ = numberIterations_;
5236	if (numberIterations_ > lastBadIteration_ + 100)
5237	moreSpecialOptions_ &= ~16; // clear check accuracy flag
5238	if (problemStatus_ < 0) {
5239	sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
5240	if (sumDualInfeasibilities_)
5241	numberDualInfeasibilities_ = 1;
5242	}
5243	#ifdef CLP_REPORT_PROGRESS
5244	if (ixxxxxx > ixxyyyy - 3) {
5245	printf("objectiveValue_ %g\n", objectiveValue_);
5246	handler_->setLogLevel(63);
5247	int nTotal = numberColumns_ + numberRows_;
5248	double newObj = 0.0;
5249	for (int i = 0; i < nTotal; i++) {
5250	if (solution_[i])
5251	newObj += solution_[i] * cost_[i];
5252	}
5253	printf("xxx obj %g\n", newObj);
5254	// for now - recompute all
5255	gutsOfSolution(NULL, NULL);
5256	newObj = 0.0;
5257	for (int i =<