Context Navigation

source: trunk/Clp/src/ClpSimplexPrimal.cpp @ 1551

Last change on this file since 1551 was 1551, checked in by mjs, 11 years ago
Rerun formatting after JJF patches (1533 and 1535).
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 161.5 KB

Line
1	/* $Id: ClpSimplexPrimal.cpp 1551 2010-05-24 23:34:44Z mjs $ */
2	// Copyright (C) 2002, International Business Machines
3	// Corporation and others. All Rights Reserved.
4
5	/* Notes on implementation of primal simplex algorithm.
6
7	When primal feasible(A):
8
9	If dual feasible, we are optimal. Otherwise choose an infeasible
10	basic variable to enter basis from a bound (B). We now need to find an
11	outgoing variable which will leave problem primal feasible so we get
12	the column of the tableau corresponding to the incoming variable
13	(with the correct sign depending if variable will go up or down).
14
15	We now perform a ratio test to determine which outgoing variable will
16	preserve primal feasibility (C). If no variable found then problem
17	is unbounded (in primal sense). If there is a variable, we then
18	perform pivot and repeat. Trivial?
19
20	-------------------------------------------
21
22	A) How do we get primal feasible? All variables have fake costs
23	outside their feasible region so it is trivial to declare problem
24	feasible. OSL did not have a phase 1/phase 2 approach but
25	instead effectively put an extra cost on infeasible basic variables
26	I am taking the same approach here, although it is generalized
27	to allow for non-linear costs and dual information.
28
29	In OSL, this weight was changed heuristically, here at present
30	it is only increased if problem looks finished. If problem is
31	feasible I check for unboundedness. If not unbounded we
32	could play with going into dual. As long as weights increase
33	any algorithm would be finite.
34
35	B) Which incoming variable to choose is a virtual base class.
36	For difficult problems steepest edge is preferred while for
37	very easy (large) problems we will need partial scan.
38
39	C) Sounds easy, but this is hardest part of algorithm.
40	1) Instead of stopping at first choice, we may be able
41	to allow that variable to go through bound and if objective
42	still improving choose again. These mini iterations can
43	increase speed by orders of magnitude but we may need to
44	go to more of a bucket choice of variable rather than looking
45	at them one by one (for speed).
46	2) Accuracy. Basic infeasibilities may be less than
47	tolerance. Pivoting on these makes objective go backwards.
48	OSL modified cost so a zero move was made, Gill et al
49	modified so a strictly positive move was made.
50	The two problems are that re-factorizations can
51	change rinfeasibilities above and below tolerances and that when
52	finished we need to reset costs and try again.
53	3) Degeneracy. Gill et al helps but may not be enough. We
54	may need more. Also it can improve speed a lot if we perturb
55	the rhs and bounds significantly.
56
57	References:
58	Forrest and Goldfarb, Steepest-edge simplex algorithms for
59	linear programming - Mathematical Programming 1992
60	Forrest and Tomlin, Implementing the simplex method for
61	the Optimization Subroutine Library - IBM Systems Journal 1992
62	Gill, Murray, Saunders, Wright A Practical Anti-Cycling
63	Procedure for Linear and Nonlinear Programming SOL report 1988
64
65
66	TODO:
67
68	a) Better recovery procedures. At present I never check on forward
69	progress. There is checkpoint/restart with reducing
70	re-factorization frequency, but this is only on singular
71	factorizations.
72	b) Fast methods for large easy problems (and also the option for
73	the code to automatically choose which method).
74	c) We need to be able to stop in various ways for OSI - this
75	is fairly easy.
76
77	*/
78
79
80	#include "CoinPragma.hpp"
81
82	#include <math.h>
83
84	#include "CoinHelperFunctions.hpp"
85	#include "ClpSimplexPrimal.hpp"
86	#include "ClpFactorization.hpp"
87	#include "ClpNonLinearCost.hpp"
88	#include "CoinPackedMatrix.hpp"
89	#include "CoinIndexedVector.hpp"
90	#include "ClpPrimalColumnPivot.hpp"
91	#include "ClpMessage.hpp"
92	#include "ClpEventHandler.hpp"
93	#include <cfloat>
94	#include <cassert>
95	#include <string>
96	#include <stdio.h>
97	#include <iostream>
98	// primal
99	int ClpSimplexPrimal::primal (int ifValuesPass , int startFinishOptions)
100	{
101
102	/*
103	Method
104
105	It tries to be a single phase approach with a weight of 1.0 being
106	given to getting optimal and a weight of infeasibilityCost_ being
107	given to getting primal feasible. In this version I have tried to
108	be clever in a stupid way. The idea of fake bounds in dual
109	seems to work so the primal analogue would be that of getting
110	bounds on reduced costs (by a presolve approach) and using
111	these for being above or below feasible region. I decided to waste
112	memory and keep these explicitly. This allows for non-linear
113	costs!
114
115	The code is designed to take advantage of sparsity so arrays are
116	seldom zeroed out from scratch or gone over in their entirety.
117	The only exception is a full scan to find incoming variable for
118	Dantzig row choice. For steepest edge we keep an updated list
119	of dual infeasibilities (actually squares).
120	On easy problems we don't need full scan - just
121	pick first reasonable.
122
123	One problem is how to tackle degeneracy and accuracy. At present
124	I am using the modification of costs which I put in OSL and which was
125	extended by Gill et al. I am still not sure of the exact details.
126
127	The flow of primal is three while loops as follows:
128
129	while (not finished) {
130
131	while (not clean solution) {
132
133	Factorize and/or clean up solution by changing bounds so
134	primal feasible. If looks finished check fake primal bounds.
135	Repeat until status is iterating (-1) or finished (0,1,2)
136
137	}
138
139	while (status==-1) {
140
141	Iterate until no pivot in or out or time to re-factorize.
142
143	Flow is:
144
145	choose pivot column (incoming variable). if none then
146	we are primal feasible so looks as if done but we need to
147	break and check bounds etc.
148
149	Get pivot column in tableau
150
151	Choose outgoing row. If we don't find one then we look
152	primal unbounded so break and check bounds etc. (Also the
153	pivot tolerance is larger after any iterations so that may be
154	reason)
155
156	If we do find outgoing row, we may have to adjust costs to
157	keep going forwards (anti-degeneracy). Check pivot will be stable
158	and if unstable throw away iteration and break to re-factorize.
159	If minor error re-factorize after iteration.
160
161	Update everything (this may involve changing bounds on
162	variables to stay primal feasible.
163
164	}
165
166	}
167
168	At present we never check we are going forwards. I overdid that in
169	OSL so will try and make a last resort.
170
171	Needs partial scan pivot in option.
172
173	May need other anti-degeneracy measures, especially if we try and use
174	loose tolerances as a way to solve in fewer iterations.
175
176	I like idea of dynamic scaling. This gives opportunity to decouple
177	different implications of scaling for accuracy, iteration count and
178	feasibility tolerance.
179
180	*/
181
182	algorithm_ = +1;
183	moreSpecialOptions_ &= ~16; // clear check replaceColumn accuracy
184
185	// save data
186	ClpDataSave data = saveData();
187	matrix_->refresh(this); // make sure matrix okay
188
189	// Save so can see if doing after dual
190	int initialStatus = problemStatus_;
191	int initialIterations = numberIterations_;
192	int initialNegDjs = -1;
193	// initialize - maybe values pass and algorithm_ is +1
194	#if 0
195	// if so - put in any superbasic costed slacks
196	if (ifValuesPass && specialOptions_ < 0x01000000) {
197	// Get column copy
198	const CoinPackedMatrix * columnCopy = matrix();
199	const int * row = columnCopy->getIndices();
200	const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
201	const int * columnLength = columnCopy->getVectorLengths();
202	//const double * element = columnCopy->getElements();
203	int n = 0;
204	for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
205	if (columnLength[iColumn] == 1) {
206	Status status = getColumnStatus(iColumn);
207	if (status != basic && status != isFree) {
208	double value = columnActivity_[iColumn];
209	if (fabs(value - columnLower_[iColumn]) > primalTolerance_ &&
210	fabs(value - columnUpper_[iColumn]) > primalTolerance_) {
211	int iRow = row[columnStart[iColumn]];
212	if (getRowStatus(iRow) == basic) {
213	setRowStatus(iRow, superBasic);
214	setColumnStatus(iColumn, basic);
215	n++;
216	}
217	}
218	}
219	}
220	}
221	printf("%d costed slacks put in basis\n", n);
222	}
223	#endif
224	// Start can skip some things in transposeTimes
225	specialOptions_ \|= 131072;
226	if (!startup(ifValuesPass, startFinishOptions)) {
227
228	// Set average theta
229	nonLinearCost_->setAverageTheta(1.0e3);
230	int lastCleaned = 0; // last time objective or bounds cleaned up
231
232	// Say no pivot has occurred (for steepest edge and updates)
233	pivotRow_ = -2;
234
235	// This says whether to restore things etc
236	int factorType = 0;
237	if (problemStatus_ < 0 && perturbation_ < 100 && !ifValuesPass) {
238	perturb(0);
239	// Can't get here if values pass
240	assert (!ifValuesPass);
241	gutsOfSolution(NULL, NULL);
242	if (handler_->logLevel() > 2) {
243	handler_->message(CLP_SIMPLEX_STATUS, messages_)
244	<< numberIterations_ << objectiveValue();
245	handler_->printing(sumPrimalInfeasibilities_ > 0.0)
246	<< sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
247	handler_->printing(sumDualInfeasibilities_ > 0.0)
248	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
249	handler_->printing(numberDualInfeasibilitiesWithoutFree_
250	< numberDualInfeasibilities_)
251	<< numberDualInfeasibilitiesWithoutFree_;
252	handler_->message() << CoinMessageEol;
253	}
254	}
255	ClpSimplex * saveModel = NULL;
256	int stopSprint = -1;
257	int sprintPass = 0;
258	int reasonableSprintIteration = 0;
259	int lastSprintIteration = 0;
260	double lastObjectiveValue = COIN_DBL_MAX;
261	// Start check for cycles
262	progress_.fillFromModel(this);
263	progress_.startCheck();
264	/*
265	Status of problem:
266	0 - optimal
267	1 - infeasible
268	2 - unbounded
269	-1 - iterating
270	-2 - factorization wanted
271	-3 - redo checking without factorization
272	-4 - looks infeasible
273	-5 - looks unbounded
274	*/
275	while (problemStatus_ < 0) {
276	int iRow, iColumn;
277	// clear
278	for (iRow = 0; iRow < 4; iRow++) {
279	rowArray_[iRow]->clear();
280	}
281
282	for (iColumn = 0; iColumn < 2; iColumn++) {
283	columnArray_[iColumn]->clear();
284	}
285
286	// give matrix (and model costs and bounds a chance to be
287	// refreshed (normally null)
288	matrix_->refresh(this);
289	// If getting nowhere - why not give it a kick
290	#if 1
291	if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (specialOptions_ & 4) == 0
292	&& initialStatus != 10) {
293	perturb(1);
294	matrix_->rhsOffset(this, true, false);
295	}
296	#endif
297	// If we have done no iterations - special
298	if (lastGoodIteration_ == numberIterations_ && factorType)
299	factorType = 3;
300	if (saveModel) {
301	// Doing sprint
302	if (sequenceIn_ < 0 \|\| numberIterations_ >= stopSprint) {
303	problemStatus_ = -1;
304	originalModel(saveModel);
305	saveModel = NULL;
306	if (sequenceIn_ < 0 && numberIterations_ < reasonableSprintIteration &&
307	sprintPass > 100)
308	primalColumnPivot_->switchOffSprint();
309	//lastSprintIteration=numberIterations_;
310	printf("End small model\n");
311	}
312	}
313
314	// may factorize, checks if problem finished
315	statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
316	if (initialStatus == 10) {
317	// cleanup phase
318	if(initialIterations != numberIterations_) {
319	if (numberDualInfeasibilities_ > 10000 && numberDualInfeasibilities_ > 10 * initialNegDjs) {
320	// getting worse - try perturbing
321	if (perturbation_ < 101 && (specialOptions_ & 4) == 0) {
322	perturb(1);
323	matrix_->rhsOffset(this, true, false);
324	statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
325	}
326	}
327	} else {
328	// save number of negative djs
329	if (!numberPrimalInfeasibilities_)
330	initialNegDjs = numberDualInfeasibilities_;
331	// make sure weight won't be changed
332	if (infeasibilityCost_ == 1.0e10)
333	infeasibilityCost_ = 1.000001e10;
334	}
335	}
336	// See if sprint says redo because of problems
337	if (numberDualInfeasibilities_ == -776) {
338	// Need new set of variables
339	problemStatus_ = -1;
340	originalModel(saveModel);
341	saveModel = NULL;
342	//lastSprintIteration=numberIterations_;
343	printf("End small model after\n");
344	statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
345	}
346	int numberSprintIterations = 0;
347	int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
348	if (problemStatus_ == 777) {
349	// problems so do one pass with normal
350	problemStatus_ = -1;
351	originalModel(saveModel);
352	saveModel = NULL;
353	// Skip factorization
354	//statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
355	statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
356	} else if (problemStatus_ < 0 && !saveModel && numberSprintColumns && firstFree_ < 0) {
357	int numberSort = 0;
358	int numberFixed = 0;
359	int numberBasic = 0;
360	reasonableSprintIteration = numberIterations_ + 100;
361	int * whichColumns = new int[numberColumns_];
362	double * weight = new double[numberColumns_];
363	int numberNegative = 0;
364	double sumNegative = 0.0;
365	// now massage weight so all basic in plus good djs
366	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
367	double dj = dj_[iColumn];
368	switch(getColumnStatus(iColumn)) {
369
370	case basic:
371	dj = -1.0e50;
372	numberBasic++;
373	break;
374	case atUpperBound:
375	dj = -dj;
376	break;
377	case isFixed:
378	dj = 1.0e50;
379	numberFixed++;
380	break;
381	case atLowerBound:
382	dj = dj;
383	break;
384	case isFree:
385	dj = -100.0 * fabs(dj);
386	break;
387	case superBasic:
388	dj = -100.0 * fabs(dj);
389	break;
390	}
391	if (dj < -dualTolerance_ && dj > -1.0e50) {
392	numberNegative++;
393	sumNegative -= dj;
394	}
395	weight[iColumn] = dj;
396	whichColumns[iColumn] = iColumn;
397	}
398	handler_->message(CLP_SPRINT, messages_)
399	<< sprintPass << numberIterations_ - lastSprintIteration << objectiveValue() << sumNegative
400	<< numberNegative
401	<< CoinMessageEol;
402	sprintPass++;
403	lastSprintIteration = numberIterations_;
404	if (objectiveValue()*optimizationDirection_ > lastObjectiveValue - 1.0e-7 && sprintPass > 5) {
405	// switch off
406	printf("Switching off sprint\n");
407	primalColumnPivot_->switchOffSprint();
408	} else {
409	lastObjectiveValue = objectiveValue() * optimizationDirection_;
410	// sort
411	CoinSort_2(weight, weight + numberColumns_, whichColumns);
412	numberSort = CoinMin(numberColumns_ - numberFixed, numberBasic + numberSprintColumns);
413	// Sort to make consistent ?
414	std::sort(whichColumns, whichColumns + numberSort);
415	saveModel = new ClpSimplex(this, numberSort, whichColumns);
416	delete [] whichColumns;
417	delete [] weight;
418	// Skip factorization
419	//statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
420	//statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
421	stopSprint = numberIterations_ + numberSprintIterations;
422	printf("Sprint with %d columns for %d iterations\n",
423	numberSprintColumns, numberSprintIterations);
424	}
425	}
426
427	// Say good factorization
428	factorType = 1;
429
430	// Say no pivot has occurred (for steepest edge and updates)
431	pivotRow_ = -2;
432
433	// exit if victory declared
434	if (problemStatus_ >= 0)
435	break;
436
437	// test for maximum iterations
438	if (hitMaximumIterations() \|\| (ifValuesPass == 2 && firstFree_ < 0)) {
439	problemStatus_ = 3;
440	break;
441	}
442
443	if (firstFree_ < 0) {
444	if (ifValuesPass) {
445	// end of values pass
446	ifValuesPass = 0;
447	int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
448	if (status >= 0) {
449	problemStatus_ = 5;
450	secondaryStatus_ = ClpEventHandler::endOfValuesPass;
451	break;
452	}
453	//#define FEB_TRY
454	#ifdef FEB_TRY
455	if (perturbation_ < 100)
456	perturb(0);
457	#endif
458	}
459	}
460	// Check event
461	{
462	int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
463	if (status >= 0) {
464	problemStatus_ = 5;
465	secondaryStatus_ = ClpEventHandler::endOfFactorization;
466	break;
467	}
468	}
469	// Iterate
470	whileIterating(ifValuesPass ? 1 : 0);
471	}
472	}
473	// if infeasible get real values
474	//printf("XXXXY final cost %g\n",infeasibilityCost_);
475	progress_.initialWeight_ = 0.0;
476	if (problemStatus_ == 1 && secondaryStatus_ != 6) {
477	infeasibilityCost_ = 0.0;
478	createRim(1 + 4);
479	delete nonLinearCost_;
480	nonLinearCost_ = new ClpNonLinearCost(this);
481	nonLinearCost_->checkInfeasibilities(0.0);
482	sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
483	numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
484	// and get good feasible duals
485	computeDuals(NULL);
486	}
487	// Stop can skip some things in transposeTimes
488	specialOptions_ &= ~131072;
489	// clean up
490	unflag();
491	finish(startFinishOptions);
492	restoreData(data);
493	return problemStatus_;
494	}
495	/*
496	Reasons to come out:
497	-1 iterations etc
498	-2 inaccuracy
499	-3 slight inaccuracy (and done iterations)
500	-4 end of values pass and done iterations
501	+0 looks optimal (might be infeasible - but we will investigate)
502	+2 looks unbounded
503	+3 max iterations
504	*/
505	int
506	ClpSimplexPrimal::whileIterating(int valuesOption)
507	{
508	// Say if values pass
509	int ifValuesPass = (firstFree_ >= 0) ? 1 : 0;
510	int returnCode = -1;
511	int superBasicType = 1;
512	if (valuesOption > 1)
513	superBasicType = 3;
514	// status stays at -1 while iterating, >=0 finished, -2 to invert
515	// status -3 to go to top without an invert
516	while (problemStatus_ == -1) {
517	//#define CLP_DEBUG 1
518	#ifdef CLP_DEBUG
519	{
520	int i;
521	// not [1] as has information
522	for (i = 0; i < 4; i++) {
523	if (i != 1)
524	rowArray_[i]->checkClear();
525	}
526	for (i = 0; i < 2; i++) {
527	columnArray_[i]->checkClear();
528	}
529	}
530	#endif
531	#if 0
532	{
533	int iPivot;
534	double * array = rowArray_[3]->denseVector();
535	int * index = rowArray_[3]->getIndices();
536	int i;
537	for (iPivot = 0; iPivot < numberRows_; iPivot++) {
538	int iSequence = pivotVariable_[iPivot];
539	unpackPacked(rowArray_[3], iSequence);
540	factorization_->updateColumn(rowArray_[2], rowArray_[3]);
541	int number = rowArray_[3]->getNumElements();
542	for (i = 0; i < number; i++) {
543	int iRow = index[i];
544	if (iRow == iPivot)
545	assert (fabs(array[i] - 1.0) < 1.0e-4);
546	else
547	assert (fabs(array[i]) < 1.0e-4);
548	}
549	rowArray_[3]->clear();
550	}
551	}
552	#endif
553	#if 0
554	nonLinearCost_->checkInfeasibilities(primalTolerance_);
555	printf("suminf %g number %d\n", nonLinearCost_->sumInfeasibilities(),
556	nonLinearCost_->numberInfeasibilities());
557	#endif
558	#if CLP_DEBUG>2
559	// very expensive
560	if (numberIterations_ > 0 && numberIterations_ < 100 && !ifValuesPass) {
561	handler_->setLogLevel(63);
562	double saveValue = objectiveValue_;
563	double * saveRow1 = new double[numberRows_];
564	double * saveRow2 = new double[numberRows_];
565	CoinMemcpyN(rowReducedCost_, numberRows_, saveRow1);
566	CoinMemcpyN(rowActivityWork_, numberRows_, saveRow2);
567	double * saveColumn1 = new double[numberColumns_];
568	double * saveColumn2 = new double[numberColumns_];
569	CoinMemcpyN(reducedCostWork_, numberColumns_, saveColumn1);
570	CoinMemcpyN(columnActivityWork_, numberColumns_, saveColumn2);
571	gutsOfSolution(NULL, NULL, false);
572	printf("xxx %d old obj %g, recomputed %g, sum primal inf %g\n",
573	numberIterations_,
574	saveValue, objectiveValue_, sumPrimalInfeasibilities_);
575	CoinMemcpyN(saveRow1, numberRows_, rowReducedCost_);
576	CoinMemcpyN(saveRow2, numberRows_, rowActivityWork_);
577	CoinMemcpyN(saveColumn1, numberColumns_, reducedCostWork_);
578	CoinMemcpyN(saveColumn2, numberColumns_, columnActivityWork_);
579	delete [] saveRow1;
580	delete [] saveRow2;
581	delete [] saveColumn1;
582	delete [] saveColumn2;
583	objectiveValue_ = saveValue;
584	}
585	#endif
586	if (!ifValuesPass) {
587	// choose column to come in
588	// can use pivotRow_ to update weights
589	// pass in list of cost changes so can do row updates (rowArray_[1])
590	// NOTE rowArray_[0] is used by computeDuals which is a
591	// slow way of getting duals but might be used
592	primalColumn(rowArray_[1], rowArray_[2], rowArray_[3],
593	columnArray_[0], columnArray_[1]);
594	} else {
595	// in values pass
596	int sequenceIn = nextSuperBasic(superBasicType, columnArray_[0]);
597	if (valuesOption > 1)
598	superBasicType = 2;
599	if (sequenceIn < 0) {
600	// end of values pass - initialize weights etc
601	handler_->message(CLP_END_VALUES_PASS, messages_)
602	<< numberIterations_;
603	primalColumnPivot_->saveWeights(this, 5);
604	problemStatus_ = -2; // factorize now
605	pivotRow_ = -1; // say no weights update
606	returnCode = -4;
607	// Clean up
608	int i;
609	for (i = 0; i < numberRows_ + numberColumns_; i++) {
610	if (getColumnStatus(i) == atLowerBound \|\| getColumnStatus(i) == isFixed)
611	solution_[i] = lower_[i];
612	else if (getColumnStatus(i) == atUpperBound)
613	solution_[i] = upper_[i];
614	}
615	break;
616	} else {
617	// normal
618	sequenceIn_ = sequenceIn;
619	valueIn_ = solution_[sequenceIn_];
620	lowerIn_ = lower_[sequenceIn_];
621	upperIn_ = upper_[sequenceIn_];
622	dualIn_ = dj_[sequenceIn_];
623	}
624	}
625	pivotRow_ = -1;
626	sequenceOut_ = -1;
627	rowArray_[1]->clear();
628	if (sequenceIn_ >= 0) {
629	// we found a pivot column
630	assert (!flagged(sequenceIn_));
631	#ifdef CLP_DEBUG
632	if ((handler_->logLevel() & 32)) {
633	char x = isColumn(sequenceIn_) ? 'C' : 'R';
634	std::cout << "pivot column " <<
635	x << sequenceWithin(sequenceIn_) << std::endl;
636	}
637	#endif
638	#ifdef CLP_DEBUG
639	{
640	int checkSequence = -2077;
641	if (checkSequence >= 0 && checkSequence < numberRows_ + numberColumns_ && !ifValuesPass) {
642	rowArray_[2]->checkClear();
643	rowArray_[3]->checkClear();
644	double * array = rowArray_[3]->denseVector();
645	int * index = rowArray_[3]->getIndices();
646	unpackPacked(rowArray_[3], checkSequence);
647	factorization_->updateColumnForDebug(rowArray_[2], rowArray_[3]);
648	int number = rowArray_[3]->getNumElements();
649	double dualIn = cost_[checkSequence];
650	int i;
651	for (i = 0; i < number; i++) {
652	int iRow = index[i];
653	int iPivot = pivotVariable_[iRow];
654	double alpha = array[i];
655	dualIn -= alpha * cost_[iPivot];
656	}
657	printf("old dj for %d was %g, recomputed %g\n", checkSequence,
658	dj_[checkSequence], dualIn);
659	rowArray_[3]->clear();
660	if (numberIterations_ > 2000)
661	exit(1);
662	}
663	}
664	#endif
665	// do second half of iteration
666	returnCode = pivotResult(ifValuesPass);
667	if (returnCode < -1 && returnCode > -5) {
668	problemStatus_ = -2; //
669	} else if (returnCode == -5) {
670	if ((moreSpecialOptions_ & 16) == 0 && factorization_->pivots()) {
671	moreSpecialOptions_ \|= 16;
672	problemStatus_ = -2;
673	}
674	// otherwise something flagged - continue;
675	} else if (returnCode == 2) {
676	problemStatus_ = -5; // looks unbounded
677	} else if (returnCode == 4) {
678	problemStatus_ = -2; // looks unbounded but has iterated
679	} else if (returnCode != -1) {
680	assert(returnCode == 3);
681	if (problemStatus_ != 5)
682	problemStatus_ = 3;
683	}
684	} else {
685	// no pivot column
686	#ifdef CLP_DEBUG
687	if (handler_->logLevel() & 32)
688	printf("** no column pivot\n");
689	#endif
690	if (nonLinearCost_->numberInfeasibilities())
691	problemStatus_ = -4; // might be infeasible
692	// Force to re-factorize early next time
693	int numberPivots = factorization_->pivots();
694	forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
695	returnCode = 0;
696	break;
697	}
698	}
699	if (valuesOption > 1)
700	columnArray_[0]->setNumElements(0);
701	return returnCode;
702	}
703	/* Checks if finished. Updates status */
704	void
705	ClpSimplexPrimal::statusOfProblemInPrimal(int & lastCleaned, int type,
706	ClpSimplexProgress * progress,
707	bool doFactorization,
708	int ifValuesPass,
709	ClpSimplex * originalModel)
710	{
711	int dummy; // for use in generalExpanded
712	int saveFirstFree = firstFree_;
713	// number of pivots done
714	int numberPivots = factorization_->pivots();
715	if (type == 2) {
716	// trouble - restore solution
717	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
718	CoinMemcpyN(savedSolution_ + numberColumns_ ,
719	numberRows_, rowActivityWork_);
720	CoinMemcpyN(savedSolution_ ,
721	numberColumns_, columnActivityWork_);
722	// restore extra stuff
723	matrix_->generalExpanded(this, 6, dummy);
724	forceFactorization_ = 1; // a bit drastic but ..
725	pivotRow_ = -1; // say no weights update
726	changeMade_++; // say change made
727	}
728	int saveThreshold = factorization_->sparseThreshold();
729	int tentativeStatus = problemStatus_;
730	int numberThrownOut = 1; // to loop round on bad factorization in values pass
731	double lastSumInfeasibility = COIN_DBL_MAX;
732	if (numberIterations_)
733	lastSumInfeasibility = nonLinearCost_->sumInfeasibilities();
734	int nPass = 0;
735	while (numberThrownOut) {
736	int nSlackBasic = 0;
737	if (nPass) {
738	for (int i = 0; i < numberRows_; i++) {
739	if (getRowStatus(i) == basic)
740	nSlackBasic++;
741	}
742	}
743	nPass++;
744	if (problemStatus_ > -3 \|\| problemStatus_ == -4) {
745	// factorize
746	// later on we will need to recover from singularities
747	// also we could skip if first time
748	// do weights
749	// This may save pivotRow_ for use
750	if (doFactorization)
751	primalColumnPivot_->saveWeights(this, 1);
752
753	if ((type && doFactorization) \|\| nSlackBasic == numberRows_) {
754	// is factorization okay?
755	int factorStatus = internalFactorize(1);
756	if (factorStatus) {
757	if (solveType_ == 2 + 8) {
758	// say odd
759	problemStatus_ = 5;
760	return;
761	}
762	if (type != 1 \|\| largestPrimalError_ > 1.0e3
763	\|\| largestDualError_ > 1.0e3) {
764	// switch off dense
765	int saveDense = factorization_->denseThreshold();
766	factorization_->setDenseThreshold(0);
767	// Go to safe
768	factorization_->pivotTolerance(0.99);
769	// make sure will do safe factorization
770	pivotVariable_[0] = -1;
771	internalFactorize(2);
772	factorization_->setDenseThreshold(saveDense);
773	// restore extra stuff
774	matrix_->generalExpanded(this, 6, dummy);
775	} else {
776	// no - restore previous basis
777	// Keep any flagged variables
778	int i;
779	for (i = 0; i < numberRows_ + numberColumns_; i++) {
780	if (flagged(i))
781	saveStatus_[i] \|= 64; //say flagged
782	}
783	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
784	if (numberPivots <= 1) {
785	// throw out something
786	if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
787	setFlagged(sequenceIn_);
788	} else if (sequenceOut_ >= 0 && getStatus(sequenceOut_) != basic) {
789	setFlagged(sequenceOut_);
790	}
791	double newTolerance = CoinMax(0.5 + 0.499 * randomNumberGenerator_.randomDouble(), factorization_->pivotTolerance());
792	factorization_->pivotTolerance(newTolerance);
793	} else {
794	// Go to safe
795	factorization_->pivotTolerance(0.99);
796	}
797	CoinMemcpyN(savedSolution_ + numberColumns_ ,
798	numberRows_, rowActivityWork_);
799	CoinMemcpyN(savedSolution_ ,
800	numberColumns_, columnActivityWork_);
801	// restore extra stuff
802	matrix_->generalExpanded(this, 6, dummy);
803	matrix_->generalExpanded(this, 5, dummy);
804	forceFactorization_ = 1; // a bit drastic but ..
805	type = 2;
806	if (internalFactorize(2) != 0) {
807	largestPrimalError_ = 1.0e4; // force other type
808	}
809	}
810	changeMade_++; // say change made
811	}
812	}
813	if (problemStatus_ != -4)
814	problemStatus_ = -3;
815	}
816	// at this stage status is -3 or -5 if looks unbounded
817	// get primal and dual solutions
818	// put back original costs and then check
819	// createRim(4); // costs do not change
820	// May need to do more if column generation
821	dummy = 4;
822	matrix_->generalExpanded(this, 9, dummy);
823	#ifndef CLP_CAUTION
824	#define CLP_CAUTION 1
825	#endif
826	#if CLP_CAUTION
827	double lastAverageInfeasibility = sumDualInfeasibilities_ /
828	static_cast<double>(numberDualInfeasibilities_ + 10);
829	#endif
830	numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
831	double sumInfeasibility = nonLinearCost_->sumInfeasibilities();
832	int reason2 = 0;
833	#if CLP_CAUTION
834	#if CLP_CAUTION==2
835	double test2 = 1.0e5;
836	#else
837	double test2 = 1.0e-1;
838	#endif
839	if (!lastSumInfeasibility && sumInfeasibility &&
840	lastAverageInfeasibility < test2 && numberPivots > 10)
841	reason2 = 3;
842	if (lastSumInfeasibility < 1.0e-6 && sumInfeasibility > 1.0e-3 &&
843	numberPivots > 10)
844	reason2 = 4;
845	#endif
846	if (numberThrownOut)
847	reason2 = 1;
848	if ((sumInfeasibility > 1.0e7 && sumInfeasibility > 100.0 * lastSumInfeasibility
849	&& factorization_->pivotTolerance() < 0.11) \|\|
850	(largestPrimalError_ > 1.0e10 && largestDualError_ > 1.0e10))
851	reason2 = 2;
852	if (reason2) {
853	problemStatus_ = tentativeStatus;
854	doFactorization = true;
855	if (numberPivots) {
856	// go back
857	// trouble - restore solution
858	CoinMemcpyN(saveStatus_, numberColumns_ + numberRows_, status_);
859	CoinMemcpyN(savedSolution_ + numberColumns_ ,
860	numberRows_, rowActivityWork_);
861	CoinMemcpyN(savedSolution_ ,
862	numberColumns_, columnActivityWork_);
863	// restore extra stuff
864	matrix_->generalExpanded(this, 6, dummy);
865	if (reason2 < 3) {
866	// Go to safe
867	factorization_->pivotTolerance(CoinMin(0.99, 1.01 * factorization_->pivotTolerance()));
868	forceFactorization_ = 1; // a bit drastic but ..
869	} else if (forceFactorization_ < 0) {
870	forceFactorization_ = CoinMin(numberPivots / 2, 100);
871	} else {
872	forceFactorization_ = CoinMin(forceFactorization_,
873	CoinMax(3, numberPivots / 2));
874	}
875	pivotRow_ = -1; // say no weights update
876	changeMade_++; // say change made
877	if (numberPivots == 1) {
878	// throw out something
879	if (sequenceIn_ >= 0 && getStatus(sequenceIn_) != basic) {
880	setFlagged(sequenceIn_);
881	} else if (sequenceOut_ >= 0 && getStatus(sequenceOut_) != basic) {
882	setFlagged(sequenceOut_);
883	}
884	}
885	type = 2; // so will restore weights
886	if (internalFactorize(2) != 0) {
887	largestPrimalError_ = 1.0e4; // force other type
888	}
889	numberPivots = 0;
890	numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= 0));
891	assert (!numberThrownOut);
892	sumInfeasibility = nonLinearCost_->sumInfeasibilities();
893	}
894	}
895	}
896	// Double check reduced costs if no action
897	if (progress->lastIterationNumber(0) == numberIterations_) {
898	if (primalColumnPivot_->looksOptimal()) {
899	numberDualInfeasibilities_ = 0;
900	sumDualInfeasibilities_ = 0.0;
901	}
902	}
903	// If in primal and small dj give up
904	if ((specialOptions_ & 1024) != 0 && !numberPrimalInfeasibilities_ && numberDualInfeasibilities_) {
905	double average = sumDualInfeasibilities_ / (static_cast<double> (numberDualInfeasibilities_));
906	if (numberIterations_ > 300 && average < 1.0e-4) {
907	numberDualInfeasibilities_ = 0;
908	sumDualInfeasibilities_ = 0.0;
909	}
910	}
911	// Check if looping
912	int loop;
913	if (type != 2 && !ifValuesPass)
914	loop = progress->looping();
915	else
916	loop = -1;
917	if (loop >= 0) {
918	if (!problemStatus_) {
919	// declaring victory
920	numberPrimalInfeasibilities_ = 0;
921	sumPrimalInfeasibilities_ = 0.0;
922	} else {
923	problemStatus_ = loop; //exit if in loop
924	problemStatus_ = 10; // instead - try other algorithm
925	numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
926	}
927	problemStatus_ = 10; // instead - try other algorithm
928	return ;
929	} else if (loop < -1) {
930	// Is it time for drastic measures
931	if (nonLinearCost_->numberInfeasibilities() && progress->badTimes() > 5 &&
932	progress->oddState() < 10 && progress->oddState() >= 0) {
933	progress->newOddState();
934	nonLinearCost_->zapCosts();
935	}
936	// something may have changed
937	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
938	}
939	// If progress then reset costs
940	if (loop == -1 && !nonLinearCost_->numberInfeasibilities() && progress->oddState() < 0) {
941	createRim(4, false); // costs back
942	delete nonLinearCost_;
943	nonLinearCost_ = new ClpNonLinearCost(this);
944	progress->endOddState();
945	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
946	}
947	// Flag to say whether to go to dual to clean up
948	bool goToDual = false;
949	// really for free variables in
950	//if((progressFlag_&2)!=0)
951	//problemStatus_=-1;;
952	progressFlag_ = 0; //reset progress flag
953
954	handler_->message(CLP_SIMPLEX_STATUS, messages_)
955	<< numberIterations_ << nonLinearCost_->feasibleReportCost();
956	handler_->printing(nonLinearCost_->numberInfeasibilities() > 0)
957	<< nonLinearCost_->sumInfeasibilities() << nonLinearCost_->numberInfeasibilities();
958	handler_->printing(sumDualInfeasibilities_ > 0.0)
959	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
960	handler_->printing(numberDualInfeasibilitiesWithoutFree_
961	< numberDualInfeasibilities_)
962	<< numberDualInfeasibilitiesWithoutFree_;
963	handler_->message() << CoinMessageEol;
964	if (!primalFeasible()) {
965	nonLinearCost_->checkInfeasibilities(primalTolerance_);
966	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
967	nonLinearCost_->checkInfeasibilities(primalTolerance_);
968	}
969	if (nonLinearCost_->numberInfeasibilities() > 0 && !progress->initialWeight_ && !ifValuesPass && infeasibilityCost_ == 1.0e10) {
970	// first time infeasible - start up weight computation
971	double * oldDj = dj_;
972	double * oldCost = cost_;
973	int numberRows2 = numberRows_ + numberExtraRows_;
974	int numberTotal = numberRows2 + numberColumns_;
975	dj_ = new double[numberTotal];
976	cost_ = new double[numberTotal];
977	reducedCostWork_ = dj_;
978	rowReducedCost_ = dj_ + numberColumns_;
979	objectiveWork_ = cost_;
980	rowObjectiveWork_ = cost_ + numberColumns_;
981	double direction = optimizationDirection_ * objectiveScale_;
982	const double * obj = objective();
983	memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
984	int iSequence;
985	if (columnScale_)
986	for (iSequence = 0; iSequence < numberColumns_; iSequence++)
987	cost_[iSequence] = obj[iSequence] * direction * columnScale_[iSequence];
988	else
989	for (iSequence = 0; iSequence < numberColumns_; iSequence++)
990	cost_[iSequence] = obj[iSequence] * direction;
991	computeDuals(NULL);
992	int numberSame = 0;
993	int numberDifferent = 0;
994	int numberZero = 0;
995	int numberFreeSame = 0;
996	int numberFreeDifferent = 0;
997	int numberFreeZero = 0;
998	int n = 0;
999	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
1000	if (getStatus(iSequence) != basic && !flagged(iSequence)) {
1001	// not basic
1002	double distanceUp = upper_[iSequence] - solution_[iSequence];
1003	double distanceDown = solution_[iSequence] - lower_[iSequence];
1004	double feasibleDj = dj_[iSequence];
1005	double infeasibleDj = oldDj[iSequence] - feasibleDj;
1006	double value = feasibleDj * infeasibleDj;
1007	if (distanceUp > primalTolerance_) {
1008	// Check if "free"
1009	if (distanceDown > primalTolerance_) {
1010	// free
1011	if (value > dualTolerance_) {
1012	numberFreeSame++;
1013	} else if(value < -dualTolerance_) {
1014	numberFreeDifferent++;
1015	dj_[n++] = feasibleDj / infeasibleDj;
1016	} else {
1017	numberFreeZero++;
1018	}
1019	} else {
1020	// should not be negative
1021	if (value > dualTolerance_) {
1022	numberSame++;
1023	} else if(value < -dualTolerance_) {
1024	numberDifferent++;
1025	dj_[n++] = feasibleDj / infeasibleDj;
1026	} else {
1027	numberZero++;
1028	}
1029	}
1030	} else if (distanceDown > primalTolerance_) {
1031	// should not be positive
1032	if (value > dualTolerance_) {
1033	numberSame++;
1034	} else if(value < -dualTolerance_) {
1035	numberDifferent++;
1036	dj_[n++] = feasibleDj / infeasibleDj;
1037	} else {
1038	numberZero++;
1039	}
1040	}
1041	}
1042	progress->initialWeight_ = -1.0;
1043	}
1044	//printf("XXXX %d same, %d different, %d zero, -- free %d %d %d\n",
1045	// numberSame,numberDifferent,numberZero,
1046	// numberFreeSame,numberFreeDifferent,numberFreeZero);
1047	// we want most to be same
1048	if (n) {
1049	double most = 0.95;
1050	std::sort(dj_, dj_ + n);
1051	int which = static_cast<int> ((1.0 - most) * static_cast<double> (n));
1052	double take = -dj_[which] * infeasibilityCost_;
1053	//printf("XXXXZ inf cost %g take %g (range %g %g)\n",infeasibilityCost_,take,-dj_[0]infeasibilityCost_,-dj_[n-1]infeasibilityCost_);
1054	take = -dj_[0] * infeasibilityCost_;
1055	infeasibilityCost_ = CoinMin(CoinMax(1000.0 * take, 1.0e8), 1.0000001e10);;
1056	//printf("XXXX increasing weight to %g\n",infeasibilityCost_);
1057	}
1058	delete [] dj_;
1059	delete [] cost_;
1060	dj_ = oldDj;
1061	cost_ = oldCost;
1062	reducedCostWork_ = dj_;
1063	rowReducedCost_ = dj_ + numberColumns_;
1064	objectiveWork_ = cost_;
1065	rowObjectiveWork_ = cost_ + numberColumns_;
1066	if (n)
1067	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1068	}
1069	double trueInfeasibility = nonLinearCost_->sumInfeasibilities();
1070	if (!nonLinearCost_->numberInfeasibilities() && infeasibilityCost_ == 1.0e10 && !ifValuesPass && true) {
1071	// relax if default
1072	infeasibilityCost_ = CoinMin(CoinMax(100.0 * sumDualInfeasibilities_, 1.0e8), 1.00000001e10);
1073	// reset looping criterion
1074	progress->reset();
1075	trueInfeasibility = 1.123456e10;
1076	}
1077	if (trueInfeasibility > 1.0) {
1078	// If infeasibility going up may change weights
1079	double testValue = trueInfeasibility - 1.0e-4 * (10.0 + trueInfeasibility);
1080	double lastInf = progress->lastInfeasibility(1);
1081	double lastInf3 = progress->lastInfeasibility(3);
1082	double thisObj = progress->lastObjective(0);
1083	double thisInf = progress->lastInfeasibility(0);
1084	thisObj += infeasibilityCost_ * 2.0 * thisInf;
1085	double lastObj = progress->lastObjective(1);
1086	lastObj += infeasibilityCost_ * 2.0 * lastInf;
1087	double lastObj3 = progress->lastObjective(3);
1088	lastObj3 += infeasibilityCost_ * 2.0 * lastInf3;
1089	if (lastObj < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj)) - 1.0e-7
1090	&& firstFree_ < 0) {
1091	if (handler_->logLevel() == 63)
1092	printf("lastobj %g this %g force %d ", lastObj, thisObj, forceFactorization_);
1093	int maxFactor = factorization_->maximumPivots();
1094	if (maxFactor > 10) {
1095	if (forceFactorization_ < 0)
1096	forceFactorization_ = maxFactor;
1097	forceFactorization_ = CoinMax(1, (forceFactorization_ >> 2));
1098	if (handler_->logLevel() == 63)
1099	printf("Reducing factorization frequency to %d\n", forceFactorization_);
1100	}
1101	} else if (lastObj3 < thisObj - 1.0e-5 * CoinMax(fabs(thisObj), fabs(lastObj3)) - 1.0e-7
1102	&& firstFree_ < 0) {
1103	if (handler_->logLevel() == 63)
1104	printf("lastobj3 %g this3 %g `force %d ", lastObj3, thisObj, forceFactorization_);
1105	int maxFactor = factorization_->maximumPivots();
1106	if (maxFactor > 10) {
1107	if (forceFactorization_ < 0)
1108	forceFactorization_ = maxFactor;
1109	forceFactorization_ = CoinMax(1, (forceFactorization_ * 2) / 3);
1110	if (handler_->logLevel() == 63)
1111	printf("Reducing factorization frequency to %d\n", forceFactorization_);
1112	}
1113	} else if(lastInf < testValue \|\| trueInfeasibility == 1.123456e10) {
1114	if (infeasibilityCost_ < 1.0e14) {
1115	infeasibilityCost_ *= 1.5;
1116	// reset looping criterion
1117	progress->reset();
1118	if (handler_->logLevel() == 63)
1119	printf("increasing weight to %g\n", infeasibilityCost_);
1120	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1121	}
1122	}
1123	}
1124	// we may wish to say it is optimal even if infeasible
1125	bool alwaysOptimal = (specialOptions_ & 1) != 0;
1126	// give code benefit of doubt
1127	if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
1128	sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1129	// say optimal (with these bounds etc)
1130	numberDualInfeasibilities_ = 0;
1131	sumDualInfeasibilities_ = 0.0;
1132	numberPrimalInfeasibilities_ = 0;
1133	sumPrimalInfeasibilities_ = 0.0;
1134	// But check if in sprint
1135	if (originalModel) {
1136	// Carry on and re-do
1137	numberDualInfeasibilities_ = -776;
1138	}
1139	// But if real primal infeasibilities nonzero carry on
1140	if (nonLinearCost_->numberInfeasibilities()) {
1141	// most likely to happen if infeasible
1142	double relaxedToleranceP = primalTolerance_;
1143	// we can't really trust infeasibilities if there is primal error
1144	double error = CoinMin(1.0e-2, largestPrimalError_);
1145	// allow tolerance at least slightly bigger than standard
1146	relaxedToleranceP = relaxedToleranceP + error;
1147	int ninfeas = nonLinearCost_->numberInfeasibilities();
1148	double sum = nonLinearCost_->sumInfeasibilities();
1149	double average = sum / static_cast<double> (ninfeas);
1150	#ifdef COIN_DEVELOP
1151	if (handler_->logLevel() > 0)
1152	printf("nonLinearCost says infeasible %d summing to %g\n",
1153	ninfeas, sum);
1154	#endif
1155	if (average > relaxedToleranceP) {
1156	sumOfRelaxedPrimalInfeasibilities_ = sum;
1157	numberPrimalInfeasibilities_ = ninfeas;
1158	sumPrimalInfeasibilities_ = sum;
1159	#ifdef COIN_DEVELOP
1160	bool unflagged =
1161	#endif
1162	unflag();
1163	#ifdef COIN_DEVELOP
1164	if (unflagged && handler_->logLevel() > 0)
1165	printf(" - but flagged variables\n");
1166	#endif
1167	}
1168	}
1169	}
1170	// had \|\|(type==3&&problemStatus_!=-5) -- ??? why ????
1171	if ((dualFeasible() \|\| problemStatus_ == -4) && !ifValuesPass) {
1172	// see if extra helps
1173	if (nonLinearCost_->numberInfeasibilities() &&
1174	(nonLinearCost_->sumInfeasibilities() > 1.0e-3 \|\| sumOfRelaxedPrimalInfeasibilities_)
1175	&& !alwaysOptimal) {
1176	//may need infeasiblity cost changed
1177	// we can see if we can construct a ray
1178	// make up a new objective
1179	double saveWeight = infeasibilityCost_;
1180	// save nonlinear cost as we are going to switch off costs
1181	ClpNonLinearCost * nonLinear = nonLinearCost_;
1182	// do twice to make sure Primal solution has settled
1183	// put non-basics to bounds in case tolerance moved
1184	// put back original costs
1185	createRim(4);
1186	nonLinearCost_->checkInfeasibilities(0.0);
1187	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1188
1189	infeasibilityCost_ = 1.0e100;
1190	// put back original costs
1191	createRim(4);
1192	nonLinearCost_->checkInfeasibilities(primalTolerance_);
1193	// may have fixed infeasibilities - double check
1194	if (nonLinearCost_->numberInfeasibilities() == 0) {
1195	// carry on
1196	problemStatus_ = -1;
1197	infeasibilityCost_ = saveWeight;
1198	nonLinearCost_->checkInfeasibilities(primalTolerance_);
1199	} else {
1200	nonLinearCost_ = NULL;
1201	// scale
1202	int i;
1203	for (i = 0; i < numberRows_ + numberColumns_; i++)
1204	cost_[i] *= 1.0e-95;
1205	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1206	nonLinearCost_ = nonLinear;
1207	infeasibilityCost_ = saveWeight;
1208	if ((infeasibilityCost_ >= 1.0e18 \|\|
1209	numberDualInfeasibilities_ == 0) && perturbation_ == 101) {
1210	goToDual = unPerturb(); // stop any further perturbation
1211	if (nonLinearCost_->sumInfeasibilities() > 1.0e-1)
1212	goToDual = false;
1213	nonLinearCost_->checkInfeasibilities(primalTolerance_);
1214	numberDualInfeasibilities_ = 1; // carry on
1215	problemStatus_ = -1;
1216	} else if (numberDualInfeasibilities_ == 0 && largestDualError_ > 1.0e-2 &&
1217	(moreSpecialOptions_ & 256) == 0) {
1218	goToDual = true;
1219	factorization_->pivotTolerance(CoinMax(0.9, factorization_->pivotTolerance()));
1220	}
1221	if (!goToDual) {
1222	if (infeasibilityCost_ >= 1.0e20 \|\|
1223	numberDualInfeasibilities_ == 0) {
1224	// we are infeasible - use as ray
1225	delete [] ray_;
1226	ray_ = new double [numberRows_];
1227	CoinMemcpyN(dual_, numberRows_, ray_);
1228	// and get feasible duals
1229	infeasibilityCost_ = 0.0;
1230	createRim(4);
1231	nonLinearCost_->checkInfeasibilities(primalTolerance_);
1232	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1233	// so will exit
1234	infeasibilityCost_ = 1.0e30;
1235	// reset infeasibilities
1236	sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();;
1237	numberPrimalInfeasibilities_ =
1238	nonLinearCost_->numberInfeasibilities();
1239	}
1240	if (infeasibilityCost_ < 1.0e20) {
1241	infeasibilityCost_ *= 5.0;
1242	// reset looping criterion
1243	progress->reset();
1244	changeMade_++; // say change made
1245	handler_->message(CLP_PRIMAL_WEIGHT, messages_)
1246	<< infeasibilityCost_
1247	<< CoinMessageEol;
1248	// put back original costs and then check
1249	createRim(4);
1250	nonLinearCost_->checkInfeasibilities(0.0);
1251	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1252	problemStatus_ = -1; //continue
1253	goToDual = false;
1254	} else {
1255	// say infeasible
1256	problemStatus_ = 1;
1257	}
1258	}
1259	}
1260	} else {
1261	// may be optimal
1262	if (perturbation_ == 101) {
1263	goToDual = unPerturb(); // stop any further perturbation
1264	if ((numberRows_ > 20000 \|\| numberDualInfeasibilities_) && !numberTimesOptimal_)
1265	goToDual = false; // Better to carry on a bit longer
1266	lastCleaned = -1; // carry on
1267	}
1268	bool unflagged = (unflag() != 0);
1269	if ( lastCleaned != numberIterations_ \|\| unflagged) {
1270	handler_->message(CLP_PRIMAL_OPTIMAL, messages_)
1271	<< primalTolerance_
1272	<< CoinMessageEol;
1273	if (numberTimesOptimal_ < 4) {
1274	numberTimesOptimal_++;
1275	changeMade_++; // say change made
1276	if (numberTimesOptimal_ == 1) {
1277	// better to have small tolerance even if slower
1278	factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
1279	}
1280	lastCleaned = numberIterations_;
1281	if (primalTolerance_ != dblParam_[ClpPrimalTolerance])
1282	handler_->message(CLP_PRIMAL_ORIGINAL, messages_)
1283	<< CoinMessageEol;
1284	double oldTolerance = primalTolerance_;
1285	primalTolerance_ = dblParam_[ClpPrimalTolerance];
1286	#if 0
1287	double * xcost = new double[numberRows_+numberColumns_];
1288	double * xlower = new double[numberRows_+numberColumns_];
1289	double * xupper = new double[numberRows_+numberColumns_];
1290	double * xdj = new double[numberRows_+numberColumns_];
1291	double * xsolution = new double[numberRows_+numberColumns_];
1292	CoinMemcpyN(cost_, (numberRows_ + numberColumns_), xcost);
1293	CoinMemcpyN(lower_, (numberRows_ + numberColumns_), xlower);
1294	CoinMemcpyN(upper_, (numberRows_ + numberColumns_), xupper);
1295	CoinMemcpyN(dj_, (numberRows_ + numberColumns_), xdj);
1296	CoinMemcpyN(solution_, (numberRows_ + numberColumns_), xsolution);
1297	#endif
1298	// put back original costs and then check
1299	createRim(4);
1300	nonLinearCost_->checkInfeasibilities(oldTolerance);
1301	#if 0
1302	int i;
1303	for (i = 0; i < numberRows_ + numberColumns_; i++) {
1304	if (cost_[i] != xcost[i])
1305	printf("** %d old cost %g new %g sol %g\n",
1306	i, xcost[i], cost_[i], solution_[i]);
1307	if (lower_[i] != xlower[i])
1308	printf("** %d old lower %g new %g sol %g\n",
1309	i, xlower[i], lower_[i], solution_[i]);
1310	if (upper_[i] != xupper[i])
1311	printf("** %d old upper %g new %g sol %g\n",
1312	i, xupper[i], upper_[i], solution_[i]);
1313	if (dj_[i] != xdj[i])
1314	printf("** %d old dj %g new %g sol %g\n",
1315	i, xdj[i], dj_[i], solution_[i]);
1316	if (solution_[i] != xsolution[i])
1317	printf("** %d old solution %g new %g sol %g\n",
1318	i, xsolution[i], solution_[i], solution_[i]);
1319	}
1320	delete [] xcost;
1321	delete [] xupper;
1322	delete [] xlower;
1323	delete [] xdj;
1324	delete [] xsolution;
1325	#endif
1326	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1327	if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
1328	sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
1329	// say optimal (with these bounds etc)
1330	numberDualInfeasibilities_ = 0;
1331	sumDualInfeasibilities_ = 0.0;
1332	numberPrimalInfeasibilities_ = 0;
1333	sumPrimalInfeasibilities_ = 0.0;
1334	}
1335	if (dualFeasible() && !nonLinearCost_->numberInfeasibilities() && lastCleaned >= 0)
1336	problemStatus_ = 0;
1337	else
1338	problemStatus_ = -1;
1339	} else {
1340	problemStatus_ = 0; // optimal
1341	if (lastCleaned < numberIterations_) {
1342	handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
1343	<< CoinMessageEol;
1344	}
1345	}
1346	} else {
1347	if (!alwaysOptimal \|\| !sumOfRelaxedPrimalInfeasibilities_)
1348	problemStatus_ = 0; // optimal
1349	else
1350	problemStatus_ = 1; // infeasible
1351	}
1352	}
1353	} else {
1354	// see if looks unbounded
1355	if (problemStatus_ == -5) {
1356	if (nonLinearCost_->numberInfeasibilities()) {
1357	if (infeasibilityCost_ > 1.0e18 && perturbation_ == 101) {
1358	// back off weight
1359	infeasibilityCost_ = 1.0e13;
1360	// reset looping criterion
1361	progress->reset();
1362	unPerturb(); // stop any further perturbation
1363	}
1364	//we need infeasiblity cost changed
1365	if (infeasibilityCost_ < 1.0e20) {
1366	infeasibilityCost_ *= 5.0;
1367	// reset looping criterion
1368	progress->reset();
1369	changeMade_++; // say change made
1370	handler_->message(CLP_PRIMAL_WEIGHT, messages_)
1371	<< infeasibilityCost_
1372	<< CoinMessageEol;
1373	// put back original costs and then check
1374	createRim(4);
1375	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
1376	problemStatus_ = -1; //continue
1377	} else {
1378	// say infeasible
1379	problemStatus_ = 1;
1380	// we are infeasible - use as ray
1381	delete [] ray_;
1382	ray_ = new double [numberRows_];
1383	CoinMemcpyN(dual_, numberRows_, ray_);
1384	}
1385	} else {
1386	// say unbounded
1387	problemStatus_ = 2;
1388	}
1389	} else {
1390	// carry on
1391	problemStatus_ = -1;
1392	if(type == 3 && problemStatus_ != -5) {
1393	//bool unflagged =
1394	unflag();
1395	if (sumDualInfeasibilities_ < 1.0e-3 \|\|
1396	(sumDualInfeasibilities_ / static_cast<double> (numberDualInfeasibilities_)) < 1.0e-5 \|\|
1397	progress->lastIterationNumber(0) == numberIterations_) {
1398	if (!numberPrimalInfeasibilities_) {
1399	if (numberTimesOptimal_ < 4) {
1400	numberTimesOptimal_++;
1401	changeMade_++; // say change made
1402	} else {
1403	problemStatus_ = 0;
1404	secondaryStatus_ = 5;
1405	}
1406	}
1407	}
1408	}
1409	}
1410	}
1411	if (problemStatus_ == 0) {
1412	double objVal = nonLinearCost_->feasibleCost();
1413	double tol = 1.0e-10 * CoinMax(fabs(objVal), fabs(objectiveValue_)) + 1.0e-8;
1414	if (fabs(objVal - objectiveValue_) > tol) {
1415	#ifdef COIN_DEVELOP
1416	if (handler_->logLevel() > 0)
1417	printf("nonLinearCost has feasible obj of %g, objectiveValue_ is %g\n",
1418	objVal, objectiveValue_);
1419	#endif
1420	objectiveValue_ = objVal;
1421	}
1422	}
1423	// save extra stuff
1424	matrix_->generalExpanded(this, 5, dummy);
1425	if (type == 0 \|\| type == 1) {
1426	if (type != 1 \|\| !saveStatus_) {
1427	// create save arrays
1428	delete [] saveStatus_;
1429	delete [] savedSolution_;
1430	saveStatus_ = new unsigned char [numberRows_+numberColumns_];
1431	savedSolution_ = new double [numberRows_+numberColumns_];
1432	}
1433	// save arrays
1434	CoinMemcpyN(status_, numberColumns_ + numberRows_, saveStatus_);
1435	CoinMemcpyN(rowActivityWork_,
1436	numberRows_, savedSolution_ + numberColumns_);
1437	CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_);
1438	}
1439	// see if in Cbc etc
1440	bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
1441	bool disaster = false;
1442	if (disasterArea_ && inCbcOrOther && disasterArea_->check()) {
1443	disasterArea_->saveInfo();
1444	disaster = true;
1445	}
1446	if (disaster)
1447	problemStatus_ = 3;
1448	if (problemStatus_ < 0 && !changeMade_) {
1449	problemStatus_ = 4; // unknown
1450	}
1451	lastGoodIteration_ = numberIterations_;
1452	if (numberIterations_ > lastBadIteration_ + 100)
1453	moreSpecialOptions_ &= ~16; // clear check accuracy flag
1454	if (goToDual \|\| (numberIterations_ > 1000 && largestPrimalError_ > 1.0e6
1455	&& largestDualError_ > 1.0e6)) {
1456	problemStatus_ = 10; // try dual
1457	// See if second call
1458	if ((moreSpecialOptions_ & 256) != 0) {
1459	numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
1460	sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
1461	// say infeasible
1462	if (numberPrimalInfeasibilities_)
1463	problemStatus_ = 1;
1464	}
1465	}
1466	// make sure first free monotonic
1467	if (firstFree_ >= 0 && saveFirstFree >= 0) {
1468	firstFree_ = (numberIterations_) ? saveFirstFree : -1;
1469	nextSuperBasic(1, NULL);
1470	}
1471	if (doFactorization) {
1472	// restore weights (if saved) - also recompute infeasibility list
1473	if (tentativeStatus > -3)
1474	primalColumnPivot_->saveWeights(this, (type < 2) ? 2 : 4);
1475	else
1476	primalColumnPivot_->saveWeights(this, 3);
1477	if (saveThreshold) {
1478	// use default at present
1479	factorization_->sparseThreshold(0);
1480	factorization_->goSparse();
1481	}
1482	}
1483	// Allow matrices to be sorted etc
1484	int fake = -999; // signal sort
1485	matrix_->correctSequence(this, fake, fake);
1486	}
1487	/*
1488	Row array has pivot column
1489	This chooses pivot row.
1490	For speed, we may need to go to a bucket approach when many
1491	variables go through bounds
1492	On exit rhsArray will have changes in costs of basic variables
1493	*/
1494	void
1495	ClpSimplexPrimal::primalRow(CoinIndexedVector * rowArray,
1496	CoinIndexedVector * rhsArray,
1497	CoinIndexedVector * spareArray,
1498	int valuesPass)
1499	{
1500	double saveDj = dualIn_;
1501	if (valuesPass && objective_->type() < 2) {
1502	dualIn_ = cost_[sequenceIn_];
1503
1504	double * work = rowArray->denseVector();
1505	int number = rowArray->getNumElements();
1506	int * which = rowArray->getIndices();
1507
1508	int iIndex;
1509	for (iIndex = 0; iIndex < number; iIndex++) {
1510
1511	int iRow = which[iIndex];
1512	double alpha = work[iIndex];
1513	int iPivot = pivotVariable_[iRow];
1514	dualIn_ -= alpha * cost_[iPivot];
1515	}
1516	// determine direction here
1517	if (dualIn_ < -dualTolerance_) {
1518	directionIn_ = 1;
1519	} else if (dualIn_ > dualTolerance_) {
1520	directionIn_ = -1;
1521	} else {
1522	// towards nearest bound
1523	if (valueIn_ - lowerIn_ < upperIn_ - valueIn_) {
1524	directionIn_ = -1;
1525	dualIn_ = dualTolerance_;
1526	} else {
1527	directionIn_ = 1;
1528	dualIn_ = -dualTolerance_;
1529	}
1530	}
1531	}
1532
1533	// sequence stays as row number until end
1534	pivotRow_ = -1;
1535	int numberRemaining = 0;
1536
1537	double totalThru = 0.0; // for when variables flip
1538	// Allow first few iterations to take tiny
1539	double acceptablePivot = 1.0e-1 * acceptablePivot_;
1540	if (numberIterations_ > 100)
1541	acceptablePivot = acceptablePivot_;
1542	if (factorization_->pivots() > 10)
1543	acceptablePivot = 1.0e+3 * acceptablePivot_; // if we have iterated be more strict
1544	else if (factorization_->pivots() > 5)
1545	acceptablePivot = 1.0e+2 * acceptablePivot_; // if we have iterated be slightly more strict
1546	else if (factorization_->pivots())
1547	acceptablePivot = acceptablePivot_; // relax
1548	double bestEverPivot = acceptablePivot;
1549	int lastPivotRow = -1;
1550	double lastPivot = 0.0;
1551	double lastTheta = 1.0e50;
1552
1553	// use spareArrays to put ones looked at in
1554	// First one is list of candidates
1555	// We could compress if we really know we won't need any more
1556	// Second array has current set of pivot candidates
1557	// with a backup list saved in double * part of indexed vector
1558
1559	// pivot elements
1560	double * spare;
1561	// indices
1562	int * index;
1563	spareArray->clear();
1564	spare = spareArray->denseVector();
1565	index = spareArray->getIndices();
1566
1567	// we also need somewhere for effective rhs
1568	double * rhs = rhsArray->denseVector();
1569	// and we can use indices to point to alpha
1570	// that way we can store fabs(alpha)
1571	int * indexPoint = rhsArray->getIndices();
1572	//int numberFlip=0; // Those which may change if flips
1573
1574	/*
1575	First we get a list of possible pivots. We can also see if the
1576	problem looks unbounded.
1577
1578	At first we increase theta and see what happens. We start
1579	theta at a reasonable guess. If in right area then we do bit by bit.
1580	We save possible pivot candidates
1581
1582	*/
1583
1584	// do first pass to get possibles
1585	// We can also see if unbounded
1586
1587	double * work = rowArray->denseVector();
1588	int number = rowArray->getNumElements();
1589	int * which = rowArray->getIndices();
1590
1591	// we need to swap sign if coming in from ub
1592	double way = directionIn_;
1593	double maximumMovement;
1594	if (way > 0.0)
1595	maximumMovement = CoinMin(1.0e30, upperIn_ - valueIn_);
1596	else
1597	maximumMovement = CoinMin(1.0e30, valueIn_ - lowerIn_);
1598
1599	double averageTheta = nonLinearCost_->averageTheta();
1600	double tentativeTheta = CoinMin(10.0 * averageTheta, maximumMovement);
1601	double upperTheta = maximumMovement;
1602	if (tentativeTheta > 0.5 * maximumMovement)
1603	tentativeTheta = maximumMovement;
1604	bool thetaAtMaximum = tentativeTheta == maximumMovement;
1605	// In case tiny bounds increase
1606	if (maximumMovement < 1.0)
1607	tentativeTheta *= 1.1;
1608	double dualCheck = fabs(dualIn_);
1609	// but make a bit more pessimistic
1610	dualCheck = CoinMax(dualCheck - 100.0 * dualTolerance_, 0.99 * dualCheck);
1611
1612	int iIndex;
1613	int pivotOne = -1;
1614	//#define CLP_DEBUG
1615	#ifdef CLP_DEBUG
1616	if (numberIterations_ == -3839 \|\| numberIterations_ == -3840) {
1617	double dj = cost_[sequenceIn_];
1618	printf("cost in on %d is %g, dual in %g\n", sequenceIn_, dj, dualIn_);
1619	for (iIndex = 0; iIndex < number; iIndex++) {
1620
1621	int iRow = which[iIndex];
1622	double alpha = work[iIndex];
1623	int iPivot = pivotVariable_[iRow];
1624	dj -= alpha * cost_[iPivot];
1625	printf("row %d var %d current %g %g %g, alpha %g so sol => %g (cost %g, dj %g)\n",
1626	iRow, iPivot, lower_[iPivot], solution_[iPivot], upper_[iPivot],
1627	alpha, solution_[iPivot] - 1.0e9 * alpha, cost_[iPivot], dj);
1628	}
1629	}
1630	#endif
1631	while (true) {
1632	pivotOne = -1;
1633	totalThru = 0.0;
1634	// We also re-compute reduced cost
1635	numberRemaining = 0;
1636	dualIn_ = cost_[sequenceIn_];
1637	#ifndef NDEBUG
1638	double tolerance = primalTolerance_ * 1.002;
1639	#endif
1640	for (iIndex = 0; iIndex < number; iIndex++) {
1641
1642	int iRow = which[iIndex];
1643	double alpha = work[iIndex];
1644	int iPivot = pivotVariable_[iRow];
1645	if (cost_[iPivot])
1646	dualIn_ -= alpha * cost_[iPivot];
1647	alpha *= way;
1648	double oldValue = solution_[iPivot];
1649	// get where in bound sequence
1650	// note that after this alpha is actually fabs(alpha)
1651	bool possible;
1652	// do computation same way as later on in primal
1653	if (alpha > 0.0) {
1654	// basic variable going towards lower bound
1655	double bound = lower_[iPivot];
1656	// must be exactly same as when used
1657	double change = tentativeTheta * alpha;
1658	possible = (oldValue - change) <= bound + primalTolerance_;
1659	oldValue -= bound;
1660	} else {
1661	// basic variable going towards upper bound
1662	double bound = upper_[iPivot];
1663	// must be exactly same as when used
1664	double change = tentativeTheta * alpha;
1665	possible = (oldValue - change) >= bound - primalTolerance_;
1666	oldValue = bound - oldValue;
1667	alpha = - alpha;
1668	}
1669	double value;
1670	assert (oldValue >= -tolerance);
1671	if (possible) {
1672	value = oldValue - upperTheta * alpha;
1673	if (value < -primalTolerance_ && alpha >= acceptablePivot) {
1674	upperTheta = (oldValue + primalTolerance_) / alpha;
1675	pivotOne = numberRemaining;
1676	}
1677	// add to list
1678	spare[numberRemaining] = alpha;
1679	rhs[numberRemaining] = oldValue;
1680	indexPoint[numberRemaining] = iIndex;
1681	index[numberRemaining++] = iRow;
1682	totalThru += alpha;
1683	setActive(iRow);
1684	//} else if (value<primalTolerance_*1.002) {
1685	// May change if is a flip
1686	//indexRhs[numberFlip++]=iRow;
1687	}
1688	}
1689	if (upperTheta < maximumMovement && totalThruinfeasibilityCost_ >= 1.0001 dualCheck) {
1690	// Can pivot here
1691	break;
1692	} else if (!thetaAtMaximum) {
1693	//printf("Going round with average theta of %g\n",averageTheta);
1694	tentativeTheta = maximumMovement;
1695	thetaAtMaximum = true; // seems to be odd compiler error
1696	} else {
1697	break;
1698	}
1699	}
1700	totalThru = 0.0;
1701
1702	theta_ = maximumMovement;
1703
1704	bool goBackOne = false;
1705	if (objective_->type() > 1)
1706	dualIn_ = saveDj;
1707
1708	//printf("%d remain out of %d\n",numberRemaining,number);
1709	int iTry = 0;
1710	#define MAXTRY 1000
1711	if (numberRemaining && upperTheta < maximumMovement) {
1712	// First check if previously chosen one will work
1713	if (pivotOne >= 0 && 0) {
1714	double thruCost = infeasibilityCost_ * spare[pivotOne];
1715	if (thruCost >= 0.99 * fabs(dualIn_))
1716	printf("Could pivot on %d as change %g dj %g\n",
1717	index[pivotOne], thruCost, dualIn_);
1718	double alpha = spare[pivotOne];
1719	double oldValue = rhs[pivotOne];
1720	theta_ = oldValue / alpha;
1721	pivotRow_ = pivotOne;
1722	// Stop loop
1723	iTry = MAXTRY;
1724	}
1725
1726	// first get ratio with tolerance
1727	for ( ; iTry < MAXTRY; iTry++) {
1728
1729	upperTheta = maximumMovement;
1730	int iBest = -1;
1731	for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
1732
1733	double alpha = spare[iIndex];
1734	double oldValue = rhs[iIndex];
1735	double value = oldValue - upperTheta * alpha;
1736
1737	if (value < -primalTolerance_ && alpha >= acceptablePivot) {
1738	upperTheta = (oldValue + primalTolerance_) / alpha;
1739	iBest = iIndex; // just in case weird numbers
1740	}
1741	}
1742
1743	// now look at best in this lot
1744	// But also see how infeasible small pivots will make
1745	double sumInfeasibilities = 0.0;
1746	double bestPivot = acceptablePivot;
1747	pivotRow_ = -1;
1748	for (iIndex = 0; iIndex < numberRemaining; iIndex++) {
1749
1750	int iRow = index[iIndex];
1751	double alpha = spare[iIndex];
1752	double oldValue = rhs[iIndex];
1753	double value = oldValue - upperTheta * alpha;
1754
1755	if (value <= 0 \|\| iBest == iIndex) {
1756	// how much would it cost to go thru and modify bound
1757	double trueAlpha = way * work[indexPoint[iIndex]];
1758	totalThru += nonLinearCost_->changeInCost(pivotVariable_[iRow], trueAlpha, rhs[iIndex]);
1759	setActive(iRow);
1760	if (alpha > bestPivot) {
1761	bestPivot = alpha;
1762	theta_ = oldValue / bestPivot;
1763	pivotRow_ = iIndex;
1764	} else if (alpha < acceptablePivot) {
1765	if (value < -primalTolerance_)
1766	sumInfeasibilities += -value - primalTolerance_;
1767	}
1768	}
1769	}
1770	if (bestPivot < 0.1 * bestEverPivot &&
1771	bestEverPivot > 1.0e-6 && bestPivot < 1.0e-3) {
1772	// back to previous one
1773	goBackOne = true;
1774	break;
1775	} else if (pivotRow_ == -1 && upperTheta > largeValue_) {
1776	if (lastPivot > acceptablePivot) {
1777	// back to previous one
1778	goBackOne = true;
1779	} else {
1780	// can only get here if all pivots so far too small
1781	}
1782	break;
1783	} else if (totalThru >= dualCheck) {
1784	if (sumInfeasibilities > primalTolerance_ && !nonLinearCost_->numberInfeasibilities()) {
1785	// Looks a bad choice
1786	if (lastPivot > acceptablePivot) {
1787	goBackOne = true;
1788	} else {
1789	// say no good
1790	dualIn_ = 0.0;
1791	}
1792	}
1793	break; // no point trying another loop
1794	} else {
1795	lastPivotRow = pivotRow_;
1796	lastTheta = theta_;
1797	if (bestPivot > bestEverPivot)
1798	bestEverPivot = bestPivot;
1799	}
1800	}
1801	// can get here without pivotRow_ set but with lastPivotRow
1802	if (goBackOne \|\| (pivotRow_ < 0 && lastPivotRow >= 0)) {
1803	// back to previous one
1804	pivotRow_ = lastPivotRow;
1805	theta_ = lastTheta;
1806	}
1807	} else if (pivotRow_ < 0 && maximumMovement > 1.0e20) {
1808	// looks unbounded
1809	valueOut_ = COIN_DBL_MAX; // say odd
1810	if (nonLinearCost_->numberInfeasibilities()) {
1811	// but infeasible??
1812	// move variable but don't pivot
1813	tentativeTheta = 1.0e50;
1814	for (iIndex = 0; iIndex < number; iIndex++) {
1815	int iRow = which[iIndex];
1816	double alpha = work[iIndex];
1817	int iPivot = pivotVariable_[iRow];
1818	alpha *= way;
1819	double oldValue = solution_[iPivot];
1820	// get where in bound sequence
1821	// note that after this alpha is actually fabs(alpha)
1822	if (alpha > 0.0) {
1823	// basic variable going towards lower bound
1824	double bound = lower_[iPivot];
1825	oldValue -= bound;
1826	} else {
1827	// basic variable going towards upper bound
1828	double bound = upper_[iPivot];
1829	oldValue = bound - oldValue;
1830	alpha = - alpha;
1831	}
1832	if (oldValue - tentativeTheta * alpha < 0.0) {
1833	tentativeTheta = oldValue / alpha;
1834	}
1835	}
1836	// If free in then see if we can get to 0.0
1837	if (lowerIn_ < -1.0e20 && upperIn_ > 1.0e20) {
1838	if (dualIn_ * valueIn_ > 0.0) {
1839	if (fabs(valueIn_) < 1.0e-2 && (tentativeTheta < fabs(valueIn_) \|\| tentativeTheta > 1.0e20)) {
1840	tentativeTheta = fabs(valueIn_);
1841	}
1842	}
1843	}
1844	if (tentativeTheta < 1.0e10)
1845	valueOut_ = valueIn_ + way * tentativeTheta;
1846	}
1847	}
1848	//if (iTry>50)
1849	//printf("** %d tries\n",iTry);
1850	if (pivotRow_ >= 0) {
1851	int position = pivotRow_; // position in list
1852	pivotRow_ = index[position];
1853	alpha_ = work[indexPoint[position]];
1854	// translate to sequence
1855	sequenceOut_ = pivotVariable_[pivotRow_];
1856	valueOut_ = solution(sequenceOut_);
1857	lowerOut_ = lower_[sequenceOut_];
1858	upperOut_ = upper_[sequenceOut_];
1859	#define MINIMUMTHETA 1.0e-12
1860	// Movement should be minimum for anti-degeneracy - unless
1861	// fixed variable out
1862	double minimumTheta;
1863	if (upperOut_ > lowerOut_)
1864	minimumTheta = MINIMUMTHETA;
1865	else
1866	minimumTheta = 0.0;
1867	// But can't go infeasible
1868	double distance;
1869	if (alpha_ * way > 0.0)
1870	distance = valueOut_ - lowerOut_;
1871	else
1872	distance = upperOut_ - valueOut_;
1873	if (distance - minimumTheta * fabs(alpha_) < -primalTolerance_)
1874	minimumTheta = CoinMax(0.0, (distance + 0.5 * primalTolerance_) / fabs(alpha_));
1875	// will we need to increase tolerance
1876	//#define CLP_DEBUG
1877	double largestInfeasibility = primalTolerance_;
1878	if (theta_ < minimumTheta && (specialOptions_ & 4) == 0 && !valuesPass) {
1879	theta_ = minimumTheta;
1880	for (iIndex = 0; iIndex < numberRemaining - numberRemaining; iIndex++) {
1881	largestInfeasibility = CoinMax(largestInfeasibility,
1882	-(rhs[iIndex] - spare[iIndex] * theta_));
1883	}
1884	//#define CLP_DEBUG
1885	#ifdef CLP_DEBUG
1886	if (largestInfeasibility > primalTolerance_ && (handler_->logLevel() & 32) > -1)
1887	printf("Primal tolerance increased from %g to %g\n",
1888	primalTolerance_, largestInfeasibility);
1889	#endif
1890	//#undef CLP_DEBUG
1891	primalTolerance_ = CoinMax(primalTolerance_, largestInfeasibility);
1892	}
1893	// Need to look at all in some cases
1894	if (theta_ > tentativeTheta) {
1895	for (iIndex = 0; iIndex < number; iIndex++)
1896	setActive(which[iIndex]);
1897	}
1898	if (way < 0.0)
1899	theta_ = - theta_;
1900	double newValue = valueOut_ - theta_ * alpha_;
1901	// If 4 bit set - Force outgoing variables to exact bound (primal)
1902	if (alpha_ * way < 0.0) {
1903	directionOut_ = -1; // to upper bound
1904	if (fabs(theta_) > 1.0e-6 \|\| (specialOptions_ & 4) != 0) {
1905	upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
1906	} else {
1907	upperOut_ = newValue;
1908	}
1909	} else {
1910	directionOut_ = 1; // to lower bound
1911	if (fabs(theta_) > 1.0e-6 \|\| (specialOptions_ & 4) != 0) {
1912	lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
1913	} else {
1914	lowerOut_ = newValue;
1915	}
1916	}
1917	dualOut_ = reducedCost(sequenceOut_);
1918	} else if (maximumMovement < 1.0e20) {
1919	// flip
1920	pivotRow_ = -2; // so we can tell its a flip
1921	sequenceOut_ = sequenceIn_;
1922	valueOut_ = valueIn_;
1923	dualOut_ = dualIn_;
1924	lowerOut_ = lowerIn_;
1925	upperOut_ = upperIn_;
1926	alpha_ = 0.0;
1927	if (way < 0.0) {
1928	directionOut_ = 1; // to lower bound
1929	theta_ = lowerOut_ - valueOut_;
1930	} else {
1931	directionOut_ = -1; // to upper bound
1932	theta_ = upperOut_ - valueOut_;
1933	}
1934	}
1935
1936	double theta1 = CoinMax(theta_, 1.0e-12);
1937	double theta2 = numberIterations_ * nonLinearCost_->averageTheta();
1938	// Set average theta
1939	nonLinearCost_->setAverageTheta((theta1 + theta2) / (static_cast<double> (numberIterations_ + 1)));
1940	//if (numberIterations_%1000==0)
1941	//printf("average theta is %g\n",nonLinearCost_->averageTheta());
1942
1943	// clear arrays
1944
1945	CoinZeroN(spare, numberRemaining);
1946
1947	// put back original bounds etc
1948	CoinMemcpyN(index, numberRemaining, rhsArray->getIndices());
1949	rhsArray->setNumElements(numberRemaining);
1950	rhsArray->setPacked();
1951	nonLinearCost_->goBackAll(rhsArray);
1952	rhsArray->clear();
1953
1954	}
1955	/*
1956	Chooses primal pivot column
1957	updateArray has cost updates (also use pivotRow_ from last iteration)
1958	Would be faster with separate region to scan
1959	and will have this (with square of infeasibility) when steepest
1960	For easy problems we can just choose one of the first columns we look at
1961	*/
1962	void
1963	ClpSimplexPrimal::primalColumn(CoinIndexedVector * updates,
1964	CoinIndexedVector * spareRow1,
1965	CoinIndexedVector * spareRow2,
1966	CoinIndexedVector * spareColumn1,
1967	CoinIndexedVector * spareColumn2)
1968	{
1969
1970	ClpMatrixBase * saveMatrix = matrix_;
1971	double * saveRowScale = rowScale_;
1972	if (scaledMatrix_) {
1973	rowScale_ = NULL;
1974	matrix_ = scaledMatrix_;
1975	}
1976	sequenceIn_ = primalColumnPivot_->pivotColumn(updates, spareRow1,
1977	spareRow2, spareColumn1,
1978	spareColumn2);
1979	if (scaledMatrix_) {
1980	matrix_ = saveMatrix;
1981	rowScale_ = saveRowScale;
1982	}
1983	if (sequenceIn_ >= 0) {
1984	valueIn_ = solution_[sequenceIn_];
1985	dualIn_ = dj_[sequenceIn_];
1986	if (nonLinearCost_->lookBothWays()) {
1987	// double check
1988	ClpSimplex::Status status = getStatus(sequenceIn_);
1989
1990	switch(status) {
1991	case ClpSimplex::atUpperBound:
1992	if (dualIn_ < 0.0) {
1993	// move to other side
1994	printf("For %d U (%g, %g, %g) dj changed from %g",
1995	sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
1996	upper_[sequenceIn_], dualIn_);
1997	dualIn_ -= nonLinearCost_->changeUpInCost(sequenceIn_);
1998	printf(" to %g\n", dualIn_);
1999	nonLinearCost_->setOne(sequenceIn_, upper_[sequenceIn_] + 2.0 * currentPrimalTolerance());
2000	setStatus(sequenceIn_, ClpSimplex::atLowerBound);
2001	}
2002	break;
2003	case ClpSimplex::atLowerBound:
2004	if (dualIn_ > 0.0) {
2005	// move to other side
2006	printf("For %d L (%g, %g, %g) dj changed from %g",
2007	sequenceIn_, lower_[sequenceIn_], solution_[sequenceIn_],
2008	upper_[sequenceIn_], dualIn_);
2009	dualIn_ -= nonLinearCost_->changeDownInCost(sequenceIn_);
2010	printf(" to %g\n", dualIn_);
2011	nonLinearCost_->setOne(sequenceIn_, lower_[sequenceIn_] - 2.0 * currentPrimalTolerance());
2012	setStatus(sequenceIn_, ClpSimplex::atUpperBound);
2013	}
2014	break;
2015	default:
2016	break;
2017	}
2018	}
2019	lowerIn_ = lower_[sequenceIn_];
2020	upperIn_ = upper_[sequenceIn_];
2021	if (dualIn_ > 0.0)
2022	directionIn_ = -1;
2023	else
2024	directionIn_ = 1;
2025	} else {
2026	sequenceIn_ = -1;
2027	}
2028	}
2029	/* The primals are updated by the given array.
2030	Returns number of infeasibilities.
2031	After rowArray will have list of cost changes
2032	*/
2033	int
2034	ClpSimplexPrimal::updatePrimalsInPrimal(CoinIndexedVector * rowArray,
2035	double theta,
2036	double & objectiveChange,
2037	int valuesPass)
2038	{
2039	// Cost on pivot row may change - may need to change dualIn
2040	double oldCost = 0.0;
2041	if (pivotRow_ >= 0)
2042	oldCost = cost_[sequenceOut_];
2043	//rowArray->scanAndPack();
2044	double * work = rowArray->denseVector();
2045	int number = rowArray->getNumElements();
2046	int * which = rowArray->getIndices();
2047
2048	int newNumber = 0;
2049	int pivotPosition = -1;
2050	nonLinearCost_->setChangeInCost(0.0);
2051	//printf("XX 4138 sol %g lower %g upper %g cost %g status %d\n",
2052	// solution_[4138],lower_[4138],upper_[4138],cost_[4138],status_[4138]);
2053	// allow for case where bound+tolerance == bound
2054	//double tolerance = 0.999999*primalTolerance_;
2055	double relaxedTolerance = 1.001 * primalTolerance_;
2056	int iIndex;
2057	if (!valuesPass) {
2058	for (iIndex = 0; iIndex < number; iIndex++) {
2059
2060	int iRow = which[iIndex];
2061	double alpha = work[iIndex];
2062	work[iIndex] = 0.0;
2063	int iPivot = pivotVariable_[iRow];
2064	double change = theta * alpha;
2065	double value = solution_[iPivot] - change;
2066	solution_[iPivot] = value;
2067	#ifndef NDEBUG
2068	// check if not active then okay
2069	if (!active(iRow) && (specialOptions_ & 4) == 0 && pivotRow_ != -1) {
2070	// But make sure one going out is feasible
2071	if (change > 0.0) {
2072	// going down
2073	if (value <= lower_[iPivot] + primalTolerance_) {
2074	if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
2075	value = lower_[iPivot];
2076	double difference = nonLinearCost_->setOne(iPivot, value);
2077	assert (!difference \|\| fabs(change) > 1.0e9);
2078	}
2079	} else {
2080	// going up
2081	if (value >= upper_[iPivot] - primalTolerance_) {
2082	if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
2083	value = upper_[iPivot];
2084	double difference = nonLinearCost_->setOne(iPivot, value);
2085	assert (!difference \|\| fabs(change) > 1.0e9);
2086	}
2087	}
2088	}
2089	#endif
2090	if (active(iRow) \|\| theta_ < 0.0) {
2091	clearActive(iRow);
2092	// But make sure one going out is feasible
2093	if (change > 0.0) {
2094	// going down
2095	if (value <= lower_[iPivot] + primalTolerance_) {
2096	if (iPivot == sequenceOut_ && value >= lower_[iPivot] - relaxedTolerance)
2097	value = lower_[iPivot];
2098	double difference = nonLinearCost_->setOne(iPivot, value);
2099	if (difference) {
2100	if (iRow == pivotRow_)
2101	pivotPosition = newNumber;
2102	work[newNumber] = difference;
2103	//change reduced cost on this
2104	dj_[iPivot] = -difference;
2105	which[newNumber++] = iRow;
2106	}
2107	}
2108	} else {
2109	// going up
2110	if (value >= upper_[iPivot] - primalTolerance_) {
2111	if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
2112	value = upper_[iPivot];
2113	double difference = nonLinearCost_->setOne(iPivot, value);
2114	if (difference) {
2115	if (iRow == pivotRow_)
2116	pivotPosition = newNumber;
2117	work[newNumber] = difference;
2118	//change reduced cost on this
2119	dj_[iPivot] = -difference;
2120	which[newNumber++] = iRow;
2121	}
2122	}
2123	}
2124	}
2125	}
2126	} else {
2127	// values pass so look at all
2128	for (iIndex = 0; iIndex < number; iIndex++) {
2129
2130	int iRow = which[iIndex];
2131	double alpha = work[iIndex];
2132	work[iIndex] = 0.0;
2133	int iPivot = pivotVariable_[iRow];
2134	double change = theta * alpha;
2135	double value = solution_[iPivot] - change;
2136	solution_[iPivot] = value;
2137	clearActive(iRow);
2138	// But make sure one going out is feasible
2139	if (change > 0.0) {
2140	// going down
2141	if (value <= lower_[iPivot] + primalTolerance_) {
2142	if (iPivot == sequenceOut_ && value > lower_[iPivot] - relaxedTolerance)
2143	value = lower_[iPivot];
2144	double difference = nonLinearCost_->setOne(iPivot, value);
2145	if (difference) {
2146	if (iRow == pivotRow_)
2147	pivotPosition = newNumber;
2148	work[newNumber] = difference;
2149	//change reduced cost on this
2150	dj_[iPivot] = -difference;
2151	which[newNumber++] = iRow;
2152	}
2153	}
2154	} else {
2155	// going up
2156	if (value >= upper_[iPivot] - primalTolerance_) {
2157	if (iPivot == sequenceOut_ && value < upper_[iPivot] + relaxedTolerance)
2158	value = upper_[iPivot];
2159	double difference = nonLinearCost_->setOne(iPivot, value);
2160	if (difference) {
2161	if (iRow == pivotRow_)
2162	pivotPosition = newNumber;
2163	work[newNumber] = difference;
2164	//change reduced cost on this
2165	dj_[iPivot] = -difference;
2166	which[newNumber++] = iRow;
2167	}
2168	}
2169	}
2170	}
2171	}
2172	objectiveChange += nonLinearCost_->changeInCost();
2173	rowArray->setPacked();
2174	#if 0
2175	rowArray->setNumElements(newNumber);
2176	rowArray->expand();
2177	if (pivotRow_ >= 0) {
2178	dualIn_ += (oldCost - cost_[sequenceOut_]);
2179	// update change vector to include pivot
2180	rowArray->add(pivotRow_, -dualIn_);
2181	// and convert to packed
2182	rowArray->scanAndPack();
2183	} else {
2184	// and convert to packed
2185	rowArray->scanAndPack();
2186	}
2187	#else
2188	if (pivotRow_ >= 0) {
2189	double dualIn = dualIn_ + (oldCost - cost_[sequenceOut_]);
2190	// update change vector to include pivot
2191	if (pivotPosition >= 0) {
2192	work[pivotPosition] -= dualIn;
2193	} else {
2194	work[newNumber] = -dualIn;
2195	which[newNumber++] = pivotRow_;
2196	}
2197	}
2198	rowArray->setNumElements(newNumber);
2199	#endif
2200	return 0;
2201	}
2202	// Perturbs problem
2203	void
2204	ClpSimplexPrimal::perturb(int type)
2205	{
2206	if (perturbation_ > 100)
2207	return; //perturbed already
2208	if (perturbation_ == 100)
2209	perturbation_ = 50; // treat as normal
2210	int savePerturbation = perturbation_;
2211	int i;
2212	if (!numberIterations_)
2213	cleanStatus(); // make sure status okay
2214	// Make sure feasible bounds
2215	if (nonLinearCost_)
2216	nonLinearCost_->feasibleBounds();
2217	// look at element range
2218	double smallestNegative;
2219	double largestNegative;
2220	double smallestPositive;
2221	double largestPositive;
2222	matrix_->rangeOfElements(smallestNegative, largestNegative,
2223	smallestPositive, largestPositive);
2224	smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
2225	largestPositive = CoinMax(fabs(largestNegative), largestPositive);
2226	double elementRatio = largestPositive / smallestPositive;
2227	if (!numberIterations_ && perturbation_ == 50) {
2228	// See if we need to perturb
2229	int numberTotal = CoinMax(numberRows_, numberColumns_);
2230	double * sort = new double[numberTotal];
2231	int nFixed = 0;
2232	for (i = 0; i < numberRows_; i++) {
2233	double lo = fabs(rowLower_[i]);
2234	double up = fabs(rowUpper_[i]);
2235	double value = 0.0;
2236	if (lo && lo < 1.0e20) {
2237	if (up && up < 1.0e20) {
2238	value = 0.5 * (lo + up);
2239	if (lo == up)
2240	nFixed++;
2241	} else {
2242	value = lo;
2243	}
2244	} else {
2245	if (up && up < 1.0e20)
2246	value = up;
2247	}
2248	sort[i] = value;
2249	}
2250	std::sort(sort, sort + numberRows_);
2251	int number = 1;
2252	double last = sort[0];
2253	for (i = 1; i < numberRows_; i++) {
2254	if (last != sort[i])
2255	number++;
2256	last = sort[i];
2257	}
2258	#ifdef KEEP_GOING_IF_FIXED
2259	//printf("ratio number diff rhs %g (%d %d fixed), element ratio %g\n",((double)number)/((double) numberRows_),
2260	// numberRows_,nFixed,elementRatio);
2261	#endif
2262	if (number * 4 > numberRows_ \|\| elementRatio > 1.0e12) {
2263	perturbation_ = 100;
2264	delete [] sort;
2265	return; // good enough
2266	}
2267	number = 0;
2268	#ifdef KEEP_GOING_IF_FIXED
2269	if (!integerType_) {
2270	// look at columns
2271	nFixed = 0;
2272	for (i = 0; i < numberColumns_; i++) {
2273	double lo = fabs(columnLower_[i]);
2274	double up = fabs(columnUpper_[i]);
2275	double value = 0.0;
2276	if (lo && lo < 1.0e20) {
2277	if (up && up < 1.0e20) {
2278	value = 0.5 * (lo + up);
2279	if (lo == up)
2280	nFixed++;
2281	} else {
2282	value = lo;
2283	}
2284	} else {
2285	if (up && up < 1.0e20)
2286	value = up;
2287	}
2288	sort[i] = value;
2289	}
2290	std::sort(sort, sort + numberColumns_);
2291	number = 1;
2292	last = sort[0];
2293	for (i = 1; i < numberColumns_; i++) {
2294	if (last != sort[i])
2295	number++;
2296	last = sort[i];
2297	}
2298	//printf("cratio number diff bounds %g (%d %d fixed)\n",((double)number)/((double) numberColumns_),
2299	// numberColumns_,nFixed);
2300	}
2301	#endif
2302	delete [] sort;
2303	if (number * 4 > numberColumns_) {
2304	perturbation_ = 100;
2305	return; // good enough
2306	}
2307	}
2308	// primal perturbation
2309	double perturbation = 1.0e-20;
2310	double bias = 1.0;
2311	int numberNonZero = 0;
2312	// maximum fraction of rhs/bounds to perturb
2313	double maximumFraction = 1.0e-5;
2314	if (perturbation_ >= 50) {
2315	perturbation = 1.0e-4;
2316	for (i = 0; i < numberColumns_ + numberRows_; i++) {
2317	if (upper_[i] > lower_[i] + primalTolerance_) {
2318	double lowerValue, upperValue;
2319	if (lower_[i] > -1.0e20)
2320	lowerValue = fabs(lower_[i]);
2321	else
2322	lowerValue = 0.0;
2323	if (upper_[i] < 1.0e20)
2324	upperValue = fabs(upper_[i]);
2325	else
2326	upperValue = 0.0;
2327	double value = CoinMax(fabs(lowerValue), fabs(upperValue));
2328	value = CoinMin(value, upper_[i] - lower_[i]);
2329	#if 1
2330	if (value) {
2331	perturbation += value;
2332	numberNonZero++;
2333	}
2334	#else
2335	perturbation = CoinMax(perturbation, value);
2336	#endif
2337	}
2338	}
2339	if (numberNonZero)
2340	perturbation /= static_cast<double> (numberNonZero);
2341	else
2342	perturbation = 1.0e-1;
2343	if (perturbation_ > 50 && perturbation_ < 55) {
2344	// reduce
2345	while (perturbation_ > 50) {
2346	perturbation_--;
2347	perturbation *= 0.25;
2348	bias *= 0.25;
2349	}
2350	} else if (perturbation_ >= 55 && perturbation_ < 60) {
2351	// increase
2352	while (perturbation_ > 55) {
2353	perturbation_--;
2354	perturbation *= 4.0;
2355	}
2356	perturbation_ = 50;
2357	}
2358	} else if (perturbation_ < 100) {
2359	perturbation = pow(10.0, perturbation_);
2360	// user is in charge
2361	maximumFraction = 1.0;
2362	}
2363	double largestZero = 0.0;
2364	double largest = 0.0;
2365	double largestPerCent = 0.0;
2366	bool printOut = (handler_->logLevel() == 63);
2367	printOut = false; //off
2368	// Check if all slack
2369	int number = 0;
2370	int iSequence;
2371	for (iSequence = 0; iSequence < numberRows_; iSequence++) {
2372	if (getRowStatus(iSequence) == basic)
2373	number++;
2374	}
2375	if (rhsScale_ > 100.0) {
2376	// tone down perturbation
2377	maximumFraction *= 0.1;
2378	}
2379	if (number != numberRows_)
2380	type = 1;
2381	// modify bounds
2382	// Change so at least 1.0e-5 and no more than 0.1
2383	// For now just no more than 0.1
2384	// printf("Pert type %d perturbation %g, maxF %g\n",type,perturbation,maximumFraction);
2385	// seems much slower???#define SAVE_PERT
2386	#ifdef SAVE_PERT
2387	if (2 * numberColumns_ > maximumPerturbationSize_) {
2388	delete [] perturbationArray_;
2389	maximumPerturbationSize_ = 2 * numberColumns_;
2390	perturbationArray_ = new double [maximumPerturbationSize_];
2391	for (int iColumn = 0; iColumn < maximumPerturbationSize_; iColumn++) {
2392	perturbationArray_[iColumn] = randomNumberGenerator_.randomDouble();
2393	}
2394	}
2395	#endif
2396	if (type == 1) {
2397	double tolerance = 100.0 * primalTolerance_;
2398	//double multiplier = perturbation*maximumFraction;
2399	for (iSequence = 0; iSequence < numberRows_ + numberColumns_; iSequence++) {
2400	if (getStatus(iSequence) == basic) {
2401	double lowerValue = lower_[iSequence];
2402	double upperValue = upper_[iSequence];
2403	if (upperValue > lowerValue + tolerance) {
2404	double solutionValue = solution_[iSequence];
2405	double difference = upperValue - lowerValue;
2406	difference = CoinMin(difference, perturbation);
2407	difference = CoinMin(difference, fabs(solutionValue) + 1.0);
2408	double value = maximumFraction * (difference + bias);
2409	value = CoinMin(value, 0.1);
2410	#ifndef SAVE_PERT
2411	value *= randomNumberGenerator_.randomDouble();
2412	#else
2413	value = perturbationArray_[2iSequence];
2414	#endif
2415	if (solutionValue - lowerValue <= primalTolerance_) {
2416	lower_[iSequence] -= value;
2417	} else if (upperValue - solutionValue <= primalTolerance_) {
2418	upper_[iSequence] += value;
2419	} else {
2420	#if 0
2421	if (iSequence >= numberColumns_) {
2422	// may not be at bound - but still perturb (unless free)
2423	if (upperValue > 1.0e30 && lowerValue < -1.0e30)
2424	value = 0.0;
2425	else
2426	value = - value; // as -1.0 in matrix
2427	} else {
2428	value = 0.0;
2429	}
2430	#else
2431	value = 0.0;
2432	#endif
2433	}
2434	if (value) {
2435	if (printOut)
2436	printf("col %d lower from %g to %g, upper from %g to %g\n",
2437	iSequence, lower_[iSequence], lowerValue, upper_[iSequence], upperValue);
2438	if (solutionValue) {
2439	largest = CoinMax(largest, value);
2440	if (value > (fabs(solutionValue) + 1.0)*largestPerCent)
2441	largestPerCent = value / (fabs(solutionValue) + 1.0);
2442	} else {
2443	largestZero = CoinMax(largestZero, value);
2444	}
2445	}
2446	}
2447	}
2448	}
2449	} else {
2450	double tolerance = 100.0 * primalTolerance_;
2451	for (i = 0; i < numberColumns_; i++) {
2452	double lowerValue = lower_[i], upperValue = upper_[i];
2453	if (upperValue > lowerValue + primalTolerance_) {
2454	double value = perturbation * maximumFraction;
2455	value = CoinMin(value, 0.1);
2456	#ifndef SAVE_PERT
2457	value *= randomNumberGenerator_.randomDouble();
2458	#else
2459	value = perturbationArray_[2i+1];
2460	#endif
2461	value *= randomNumberGenerator_.randomDouble();
2462	if (savePerturbation != 50) {
2463	if (fabs(value) <= primalTolerance_)
2464	value = 0.0;
2465	if (lowerValue > -1.0e20 && lowerValue)
2466	lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2467	if (upperValue < 1.0e20 && upperValue)
2468	upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2469	} else if (value) {
2470	double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2471	// get in range
2472	if (valueL <= tolerance) {
2473	valueL *= 10.0;
2474	while (valueL <= tolerance)
2475	valueL *= 10.0;
2476	} else if (valueL > 1.0) {
2477	valueL *= 0.1;
2478	while (valueL > 1.0)
2479	valueL *= 0.1;
2480	}
2481	if (lowerValue > -1.0e20 && lowerValue)
2482	lowerValue -= valueL;
2483	double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2484	// get in range
2485	if (valueU <= tolerance) {
2486	valueU *= 10.0;
2487	while (valueU <= tolerance)
2488	valueU *= 10.0;
2489	} else if (valueU > 1.0) {
2490	valueU *= 0.1;
2491	while (valueU > 1.0)
2492	valueU *= 0.1;
2493	}
2494	if (upperValue < 1.0e20 && upperValue)
2495	upperValue += valueU;
2496	}
2497	if (lowerValue != lower_[i]) {
2498	double difference = fabs(lowerValue - lower_[i]);
2499	largest = CoinMax(largest, difference);
2500	if (difference > fabs(lower_[i])*largestPerCent)
2501	largestPerCent = fabs(difference / lower_[i]);
2502	}
2503	if (upperValue != upper_[i]) {
2504	double difference = fabs(upperValue - upper_[i]);
2505	largest = CoinMax(largest, difference);
2506	if (difference > fabs(upper_[i])*largestPerCent)
2507	largestPerCent = fabs(difference / upper_[i]);
2508	}
2509	if (printOut)
2510	printf("col %d lower from %g to %g, upper from %g to %g\n",
2511	i, lower_[i], lowerValue, upper_[i], upperValue);
2512	}
2513	lower_[i] = lowerValue;
2514	upper_[i] = upperValue;
2515	}
2516	for (; i < numberColumns_ + numberRows_; i++) {
2517	double lowerValue = lower_[i], upperValue = upper_[i];
2518	double value = perturbation * maximumFraction;
2519	value = CoinMin(value, 0.1);
2520	value *= randomNumberGenerator_.randomDouble();
2521	if (upperValue > lowerValue + tolerance) {
2522	if (savePerturbation != 50) {
2523	if (fabs(value) <= primalTolerance_)
2524	value = 0.0;
2525	if (lowerValue > -1.0e20 && lowerValue)
2526	lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2527	if (upperValue < 1.0e20 && upperValue)
2528	upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2529	} else if (value) {
2530	double valueL = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2531	// get in range
2532	if (valueL <= tolerance) {
2533	valueL *= 10.0;
2534	while (valueL <= tolerance)
2535	valueL *= 10.0;
2536	} else if (valueL > 1.0) {
2537	valueL *= 0.1;
2538	while (valueL > 1.0)
2539	valueL *= 0.1;
2540	}
2541	if (lowerValue > -1.0e20 && lowerValue)
2542	lowerValue -= valueL;
2543	double valueU = value * (CoinMax(1.0e-2, 1.0e-5 * fabs(upperValue)));
2544	// get in range
2545	if (valueU <= tolerance) {
2546	valueU *= 10.0;
2547	while (valueU <= tolerance)
2548	valueU *= 10.0;
2549	} else if (valueU > 1.0) {
2550	valueU *= 0.1;
2551	while (valueU > 1.0)
2552	valueU *= 0.1;
2553	}
2554	if (upperValue < 1.0e20 && upperValue)
2555	upperValue += valueU;
2556	}
2557	} else if (upperValue > 0.0) {
2558	upperValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2559	lowerValue -= value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2560	} else if (upperValue < 0.0) {
2561	upperValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2562	lowerValue += value * (CoinMax(1.0e-2, 1.0e-5 * fabs(lowerValue)));
2563	} else {
2564	}
2565	if (lowerValue != lower_[i]) {
2566	double difference = fabs(lowerValue - lower_[i]);
2567	largest = CoinMax(largest, difference);
2568	if (difference > fabs(lower_[i])*largestPerCent)
2569	largestPerCent = fabs(difference / lower_[i]);
2570	}
2571	if (upperValue != upper_[i]) {
2572	double difference = fabs(upperValue - upper_[i]);
2573	largest = CoinMax(largest, difference);
2574	if (difference > fabs(upper_[i])*largestPerCent)
2575	largestPerCent = fabs(difference / upper_[i]);
2576	}
2577	if (printOut)
2578	printf("row %d lower from %g to %g, upper from %g to %g\n",
2579	i - numberColumns_, lower_[i], lowerValue, upper_[i], upperValue);
2580	lower_[i] = lowerValue;
2581	upper_[i] = upperValue;
2582	}
2583	}
2584	// Clean up
2585	for (i = 0; i < numberColumns_ + numberRows_; i++) {
2586	switch(getStatus(i)) {
2587
2588	case basic:
2589	break;
2590	case atUpperBound:
2591	solution_[i] = upper_[i];
2592	break;
2593	case isFixed:
2594	case atLowerBound:
2595	solution_[i] = lower_[i];
2596	break;
2597	case isFree:
2598	break;
2599	case superBasic:
2600	break;
2601	}
2602	}
2603	handler_->message(CLP_SIMPLEX_PERTURB, messages_)
2604	<< 100.0 * maximumFraction << perturbation << largest << 100.0 * largestPerCent << largestZero
2605	<< CoinMessageEol;
2606	// redo nonlinear costs
2607	// say perturbed
2608	perturbation_ = 101;
2609	}
2610	// un perturb
2611	bool
2612	ClpSimplexPrimal::unPerturb()
2613	{
2614	if (perturbation_ != 101)
2615	return false;
2616	// put back original bounds and costs
2617	createRim(1 + 4);
2618	sanityCheck();
2619	// unflag
2620	unflag();
2621	// get a valid nonlinear cost function
2622	delete nonLinearCost_;
2623	nonLinearCost_ = new ClpNonLinearCost(this);
2624	perturbation_ = 102; // stop any further perturbation
2625	// move non basic variables to new bounds
2626	nonLinearCost_->checkInfeasibilities(0.0);
2627	#if 1
2628	// Try using dual
2629	return true;
2630	#else
2631	gutsOfSolution(NULL, NULL, ifValuesPass != 0);
2632	return false;
2633	#endif
2634
2635	}
2636	// Unflag all variables and return number unflagged
2637	int
2638	ClpSimplexPrimal::unflag()
2639	{
2640	int i;
2641	int number = numberRows_ + numberColumns_;
2642	int numberFlagged = 0;
2643	// we can't really trust infeasibilities if there is dual error
2644	// allow tolerance bigger than standard to check on duals
2645	double relaxedToleranceD = dualTolerance_ + CoinMin(1.0e-2, 10.0 * largestDualError_);
2646	for (i = 0; i < number; i++) {
2647	if (flagged(i)) {
2648	clearFlagged(i);
2649	// only say if reasonable dj
2650	if (fabs(dj_[i]) > relaxedToleranceD)
2651	numberFlagged++;
2652	}
2653	}
2654	numberFlagged += matrix_->generalExpanded(this, 8, i);
2655	if (handler_->logLevel() > 2 && numberFlagged && objective_->type() > 1)
2656	printf("%d unflagged\n", numberFlagged);
2657	return numberFlagged;
2658	}
2659	// Do not change infeasibility cost and always say optimal
2660	void
2661	ClpSimplexPrimal::alwaysOptimal(bool onOff)
2662	{
2663	if (onOff)
2664	specialOptions_ \|= 1;
2665	else
2666	specialOptions_ &= ~1;
2667	}
2668	bool
2669	ClpSimplexPrimal::alwaysOptimal() const
2670	{
2671	return (specialOptions_ & 1) != 0;
2672	}
2673	// Flatten outgoing variables i.e. - always to exact bound
2674	void
2675	ClpSimplexPrimal::exactOutgoing(bool onOff)
2676	{
2677	if (onOff)
2678	specialOptions_ \|= 4;
2679	else
2680	specialOptions_ &= ~4;
2681	}
2682	bool
2683	ClpSimplexPrimal::exactOutgoing() const
2684	{
2685	return (specialOptions_ & 4) != 0;
2686	}
2687	/*
2688	Reasons to come out (normal mode/user mode):
2689	-1 normal
2690	-2 factorize now - good iteration/ NA
2691	-3 slight inaccuracy - refactorize - iteration done/ same but factor done
2692	-4 inaccuracy - refactorize - no iteration/ NA
2693	-5 something flagged - go round again/ pivot not possible
2694	+2 looks unbounded
2695	+3 max iterations (iteration done)
2696	*/
2697	int
2698	ClpSimplexPrimal::pivotResult(int ifValuesPass)
2699	{
2700
2701	bool roundAgain = true;
2702	int returnCode = -1;
2703
2704	// loop round if user setting and doing refactorization
2705	while (roundAgain) {
2706	roundAgain = false;
2707	returnCode = -1;
2708	pivotRow_ = -1;
2709	sequenceOut_ = -1;
2710	rowArray_[1]->clear();
2711	#if 0
2712	{
2713	int seq[] = {612, 643};
2714	int k;
2715	for (k = 0; k < sizeof(seq) / sizeof(int); k++) {
2716	int iSeq = seq[k];
2717	if (getColumnStatus(iSeq) != basic) {
2718	double djval;
2719	double * work;
2720	int number;
2721	int * which;
2722
2723	int iIndex;
2724	unpack(rowArray_[1], iSeq);
2725	factorization_->updateColumn(rowArray_[2], rowArray_[1]);
2726	djval = cost_[iSeq];
2727	work = rowArray_[1]->denseVector();
2728	number = rowArray_[1]->getNumElements();
2729	which = rowArray_[1]->getIndices();
2730
2731	for (iIndex = 0; iIndex < number; iIndex++) {
2732
2733	int iRow = which[iIndex];
2734	double alpha = work[iRow];
2735	int iPivot = pivotVariable_[iRow];
2736	djval -= alpha * cost_[iPivot];
2737	}
2738	double comp = 1.0e-8 + 1.0e-7 * (CoinMax(fabs(dj_[iSeq]), fabs(djval)));
2739	if (fabs(djval - dj_[iSeq]) > comp)
2740	printf("Bad dj %g for %d - true is %g\n",
2741	dj_[iSeq], iSeq, djval);
2742	assert (fabs(djval) < 1.0e-3 \|\| djval * dj_[iSeq] > 0.0);
2743	rowArray_[1]->clear();
2744	}
2745	}
2746	}
2747	#endif
2748
2749	// we found a pivot column
2750	// update the incoming column
2751	unpackPacked(rowArray_[1]);
2752	// save reduced cost
2753	double saveDj = dualIn_;
2754	factorization_->updateColumnFT(rowArray_[2], rowArray_[1]);
2755	// Get extra rows
2756	matrix_->extendUpdated(this, rowArray_[1], 0);
2757	// do ratio test and re-compute dj
2758	#ifdef CLP_USER_DRIVEN
2759	if (solveType_ != 2 \|\| (moreSpecialOptions_ & 512) == 0) {
2760	#endif
2761	primalRow(rowArray_[1], rowArray_[3], rowArray_[2],
2762	ifValuesPass);
2763	#ifdef CLP_USER_DRIVEN
2764	} else {
2765	int status = eventHandler_->event(ClpEventHandler::pivotRow);
2766	if (status >= 0) {
2767	problemStatus_ = 5;
2768	secondaryStatus_ = ClpEventHandler::pivotRow;
2769	break;
2770	}
2771	}
2772	#endif
2773	if (ifValuesPass) {
2774	saveDj = dualIn_;
2775	//assert (fabs(alpha_)>=1.0e-5\|\|(objective_->type()<2\|\|!objective_->activated())\|\|pivotRow_==-2);
2776	if (pivotRow_ == -1 \|\| (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
2777	if(fabs(dualIn_) < 1.0e2 * dualTolerance_ && objective_->type() < 2) {
2778	// try other way
2779	directionIn_ = -directionIn_;
2780	primalRow(rowArray_[1], rowArray_[3], rowArray_[2],
2781	0);
2782	}
2783	if (pivotRow_ == -1 \|\| (pivotRow_ >= 0 && fabs(alpha_) < 1.0e-5)) {
2784	if (solveType_ == 1) {
2785	// reject it
2786	char x = isColumn(sequenceIn_) ? 'C' : 'R';
2787	handler_->message(CLP_SIMPLEX_FLAG, messages_)
2788	<< x << sequenceWithin(sequenceIn_)
2789	<< CoinMessageEol;
2790	setFlagged(sequenceIn_);
2791	progress_.clearBadTimes();
2792	lastBadIteration_ = numberIterations_; // say be more cautious
2793	clearAll();
2794	pivotRow_ = -1;
2795	}
2796	returnCode = -5;
2797	break;
2798	}
2799	}
2800	}
2801	// need to clear toIndex_ in gub
2802	// ? when can I clear stuff
2803	// Clean up any gub stuff
2804	matrix_->extendUpdated(this, rowArray_[1], 1);
2805	double checkValue = 1.0e-2;
2806	if (largestDualError_ > 1.0e-5)
2807	checkValue = 1.0e-1;
2808	double test2 = dualTolerance_;
2809	double test1 = 1.0e-20;
2810	#if 0 //def FEB_TRY
2811	if (factorization_->pivots() < 1) {
2812	test1 = -1.0e-4;
2813	if ((saveDj < 0.0 && dualIn_ < -1.0e-5 * dualTolerance_) \|\|
2814	(saveDj > 0.0 && dualIn_ > 1.0e-5 * dualTolerance_))
2815	test2 = 0.0; // allow through
2816	}
2817	#endif
2818	if (!ifValuesPass && solveType_ == 1 && (saveDj * dualIn_ < test1 \|\|
2819	fabs(saveDj - dualIn_) > checkValue*(1.0 + fabs(saveDj)) \|\|
2820	fabs(dualIn_) < test2)) {
2821	if (!(saveDj * dualIn_ > 0.0 && CoinMin(fabs(saveDj), fabs(dualIn_)) >
2822	1.0e5)) {
2823	char x = isColumn(sequenceIn_) ? 'C' : 'R';
2824	handler_->message(CLP_PRIMAL_DJ, messages_)
2825	<< x << sequenceWithin(sequenceIn_)
2826	<< saveDj << dualIn_
2827	<< CoinMessageEol;
2828	if(lastGoodIteration_ != numberIterations_) {
2829	clearAll();
2830	pivotRow_ = -1; // say no weights update
2831	returnCode = -4;
2832	if(lastGoodIteration_ + 1 == numberIterations_) {
2833	// not looking wonderful - try cleaning bounds
2834	// put non-basics to bounds in case tolerance moved
2835	nonLinearCost_->checkInfeasibilities(0.0);
2836	}
2837	sequenceOut_ = -1;
2838	break;
2839	} else {
2840	// take on more relaxed criterion
2841	if (saveDj * dualIn_ < test1 \|\|
2842	fabs(saveDj - dualIn_) > 2.0e-1 * (1.0 + fabs(dualIn_)) \|\|
2843	fabs(dualIn_) < test2) {
2844	// need to reject something
2845	char x = isColumn(sequenceIn_) ? 'C' : 'R';
2846	handler_->message(CLP_SIMPLEX_FLAG, messages_)
2847	<< x << sequenceWithin(sequenceIn_)
2848	<< CoinMessageEol;
2849	setFlagged(sequenceIn_);
2850	#ifdef FEB_TRY
2851	// Make safer?
2852	factorization_->saferTolerances (1.0e-15, -1.03);
2853	#endif
2854	progress_.clearBadTimes();
2855	lastBadIteration_ = numberIterations_; // say be more cautious
2856	clearAll();
2857	pivotRow_ = -1;
2858	returnCode = -5;
2859	sequenceOut_ = -1;
2860	break;
2861	}
2862	}
2863	} else {
2864	//printf("%d %g %g\n",numberIterations_,saveDj,dualIn_);
2865	}
2866	}
2867	if (pivotRow_ >= 0) {
2868	if (solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) {
2869	// **** Coding for user interface
2870	// do ray
2871	primalRay(rowArray_[1]);
2872	// update duals
2873	// as packed need to find pivot row
2874	//assert (rowArray_[1]->packedMode());
2875	//int i;
2876
2877	//alpha_ = rowArray_[1]->denseVector()[pivotRow_];
2878	CoinAssert (fabs(alpha_) > 1.0e-8);
2879	double multiplier = dualIn_ / alpha_;
2880	rowArray_[0]->insert(pivotRow_, multiplier);
2881	factorization_->updateColumnTranspose(rowArray_[2], rowArray_[0]);
2882	// put row of tableau in rowArray[0] and columnArray[0]
2883	matrix_->transposeTimes(this, -1.0,
2884	rowArray_[0], columnArray_[1], columnArray_[0]);
2885	// update column djs
2886	int i;
2887	int * index = columnArray_[0]->getIndices();
2888	int number = columnArray_[0]->getNumElements();
2889	double * element = columnArray_[0]->denseVector();
2890	for (i = 0; i < number; i++) {
2891	int ii = index[i];
2892	dj_[ii] += element[ii];
2893	reducedCost_[ii] = dj_[ii];
2894	element[ii] = 0.0;
2895	}
2896	columnArray_[0]->setNumElements(0);
2897	// and row djs
2898	index = rowArray_[0]->getIndices();
2899	number = rowArray_[0]->getNumElements();
2900	element = rowArray_[0]->denseVector();
2901	for (i = 0; i < number; i++) {
2902	int ii = index[i];
2903	dj_[ii+numberColumns_] += element[ii];
2904	dual_[ii] = dj_[ii+numberColumns_];
2905	element[ii] = 0.0;
2906	}
2907	rowArray_[0]->setNumElements(0);
2908	// check incoming
2909	CoinAssert (fabs(dj_[sequenceIn_]) < 1.0e-1);
2910	}
2911	// if stable replace in basis
2912	// If gub or odd then alpha and pivotRow may change
2913	int updateType = 0;
2914	int updateStatus = matrix_->generalExpanded(this, 3, updateType);
2915	if (updateType >= 0)
2916	updateStatus = factorization_->replaceColumn(this,
2917	rowArray_[2],
2918	rowArray_[1],
2919	pivotRow_,
2920	alpha_,
2921	(moreSpecialOptions_ & 16) != 0);
2922
2923	// if no pivots, bad update but reasonable alpha - take and invert
2924	if (updateStatus == 2 &&
2925	lastGoodIteration_ == numberIterations_ && fabs(alpha_) > 1.0e-5)
2926	updateStatus = 4;
2927	if (updateStatus == 1 \|\| updateStatus == 4) {
2928	// slight error
2929	if (factorization_->pivots() > 5 \|\| updateStatus == 4) {
2930	returnCode = -3;
2931	}
2932	} else if (updateStatus == 2) {
2933	// major error
2934	// better to have small tolerance even if slower
2935	factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), 1.0e-15));
2936	int maxFactor = factorization_->maximumPivots();
2937	if (maxFactor > 10) {
2938	if (forceFactorization_ < 0)
2939	forceFactorization_ = maxFactor;
2940	forceFactorization_ = CoinMax(1, (forceFactorization_ >> 1));
2941	}
2942	// later we may need to unwind more e.g. fake bounds
2943	if(lastGoodIteration_ != numberIterations_) {
2944	clearAll();
2945	pivotRow_ = -1;
2946	if (solveType_ == 1 \|\| (moreSpecialOptions_ & 512) != 0) {
2947	returnCode = -4;
2948	break;
2949	} else {
2950	// user in charge - re-factorize
2951	int lastCleaned = 0;
2952	ClpSimplexProgress dummyProgress;
2953	if (saveStatus_)
2954	statusOfProblemInPrimal(lastCleaned, 1, &dummyProgress, true, ifValuesPass);
2955	else
2956	statusOfProblemInPrimal(lastCleaned, 0, &dummyProgress, true, ifValuesPass);
2957	roundAgain = true;
2958	continue;
2959	}
2960	} else {
2961	// need to reject something
2962	if (solveType_ == 1) {
2963	char x = isColumn(sequenceIn_) ? 'C' : 'R';
2964	handler_->message(CLP_SIMPLEX_FLAG, messages_)
2965	<< x << sequenceWithin(sequenceIn_)
2966	<< CoinMessageEol;
2967	setFlagged(sequenceIn_);
2968	progress_.clearBadTimes();
2969	}
2970	lastBadIteration_ = numberIterations_; // say be more cautious
2971	clearAll();
2972	pivotRow_ = -1;
2973	sequenceOut_ = -1;
2974	returnCode = -5;
2975	break;
2976
2977	}
2978	} else if (updateStatus == 3) {
2979	// out of memory
2980	// increase space if not many iterations
2981	if (factorization_->pivots() <
2982	0.5 * factorization_->maximumPivots() &&
2983	factorization_->pivots() < 200)
2984	factorization_->areaFactor(
2985	factorization_->areaFactor() * 1.1);
2986	returnCode = -2; // factorize now
2987	} else if (updateStatus == 5) {
2988	problemStatus_ = -2; // factorize now
2989	}
2990	// here do part of steepest - ready for next iteration
2991	if (!ifValuesPass)
2992	primalColumnPivot_->updateWeights(rowArray_[1]);
2993	} else {
2994	if (pivotRow_ == -1) {
2995	// no outgoing row is valid
2996	if (valueOut_ != COIN_DBL_MAX) {
2997	double objectiveChange = 0.0;
2998	theta_ = valueOut_ - valueIn_;
2999	updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, ifValuesPass);
3000	solution_[sequenceIn_] += theta_;
3001	}
3002	rowArray_[0]->clear();
3003	if (!factorization_->pivots() && acceptablePivot_ <= 1.0e-8) {
3004	returnCode = 2; //say looks unbounded
3005	// do ray
3006	primalRay(rowArray_[1]);
3007	} else if (solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) {
3008	// refactorize
3009	int lastCleaned = 0;
3010	ClpSimplexProgress dummyProgress;
3011	if (saveStatus_)
3012	statusOfProblemInPrimal(lastCleaned, 1, &dummyProgress, true, ifValuesPass);
3013	else
3014	statusOfProblemInPrimal(lastCleaned, 0, &dummyProgress, true, ifValuesPass);
3015	roundAgain = true;
3016	continue;
3017	} else {
3018	acceptablePivot_ = 1.0e-8;
3019	returnCode = 4; //say looks unbounded but has iterated
3020	}
3021	break;
3022	} else {
3023	// flipping from bound to bound
3024	}
3025	}
3026
3027	double oldCost = 0.0;
3028	if (sequenceOut_ >= 0)
3029	oldCost = cost_[sequenceOut_];
3030	// update primal solution
3031
3032	double objectiveChange = 0.0;
3033	// after this rowArray_[1] is not empty - used to update djs
3034	// If pivot row >= numberRows then may be gub
3035	int savePivot = pivotRow_;
3036	if (pivotRow_ >= numberRows_)
3037	pivotRow_ = -1;
3038	updatePrimalsInPrimal(rowArray_[1], theta_, objectiveChange, ifValuesPass);
3039	pivotRow_ = savePivot;
3040
3041	double oldValue = valueIn_;
3042	if (directionIn_ == -1) {
3043	// as if from upper bound
3044	if (sequenceIn_ != sequenceOut_) {
3045	// variable becoming basic
3046	valueIn_ -= fabs(theta_);
3047	} else {
3048	valueIn_ = lowerIn_;
3049	}
3050	} else {
3051	// as if from lower bound
3052	if (sequenceIn_ != sequenceOut_) {
3053	// variable becoming basic
3054	valueIn_ += fabs(theta_);
3055	} else {
3056	valueIn_ = upperIn_;
3057	}
3058	}
3059	objectiveChange += dualIn_ * (valueIn_ - oldValue);
3060	// outgoing
3061	if (sequenceIn_ != sequenceOut_) {
3062	if (directionOut_ > 0) {
3063	valueOut_ = lowerOut_;
3064	} else {
3065	valueOut_ = upperOut_;
3066	}
3067	if(valueOut_ < lower_[sequenceOut_] - primalTolerance_)
3068	valueOut_ = lower_[sequenceOut_] - 0.9 * primalTolerance_;
3069	else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_)
3070	valueOut_ = upper_[sequenceOut_] + 0.9 * primalTolerance_;
3071	// may not be exactly at bound and bounds may have changed
3072	// Make sure outgoing looks feasible
3073	directionOut_ = nonLinearCost_->setOneOutgoing(sequenceOut_, valueOut_);
3074	// May have got inaccurate
3075	//if (oldCost!=cost_[sequenceOut_])
3076	//printf("costchange on %d from %g to %g\n",sequenceOut_,
3077	// oldCost,cost_[sequenceOut_]);
3078	if (solveType_ != 2)
3079	dj_[sequenceOut_] = cost_[sequenceOut_] - oldCost; // normally updated next iteration
3080	solution_[sequenceOut_] = valueOut_;
3081	}
3082	// change cost and bounds on incoming if primal
3083	nonLinearCost_->setOne(sequenceIn_, valueIn_);
3084	int whatNext = housekeeping(objectiveChange);
3085	//nonLinearCost_->validate();
3086	#if CLP_DEBUG >1
3087	{
3088	double sum;
3089	int ninf = matrix_->checkFeasible(this, sum);
3090	if (ninf)
3091	printf("infeas %d\n", ninf);
3092	}
3093	#endif
3094	if (whatNext == 1) {
3095	returnCode = -2; // refactorize
3096	} else if (whatNext == 2) {
3097	// maximum iterations or equivalent
3098	returnCode = 3;
3099	} else if(numberIterations_ == lastGoodIteration_
3100	+ 2 * factorization_->maximumPivots()) {
3101	// done a lot of flips - be safe
3102	returnCode = -2; // refactorize
3103	}
3104	// Check event
3105	{
3106	int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3107	if (status >= 0) {
3108	problemStatus_ = 5;
3109	secondaryStatus_ = ClpEventHandler::endOfIteration;
3110	returnCode = 3;
3111	}
3112	}
3113	}
3114	if ((solveType_ == 2 && (moreSpecialOptions_ & 512) == 0) &&
3115	(returnCode == -2 \|\| returnCode == -3)) {
3116	// refactorize here
3117	int lastCleaned = 0;
3118	ClpSimplexProgress dummyProgress;
3119	if (saveStatus_)
3120	statusOfProblemInPrimal(lastCleaned, 1, &dummyProgress, true, ifValuesPass);
3121	else
3122	statusOfProblemInPrimal(lastCleaned, 0, &dummyProgress, true, ifValuesPass);
3123	if (problemStatus_ == 5) {
3124	printf("Singular basis\n");
3125	problemStatus_ = -1;
3126	returnCode = 5;
3127	}
3128	}
3129	#ifdef CLP_DEBUG
3130	{
3131	int i;
3132	// not [1] as may have information
3133	for (i = 0; i < 4; i++) {
3134	if (i != 1)
3135	rowArray_[i]->checkClear();
3136	}
3137	for (i = 0; i < 2; i++) {
3138	columnArray_[i]->checkClear();
3139	}
3140	}
3141	#endif
3142	return returnCode;
3143	}
3144	// Create primal ray
3145	void
3146	ClpSimplexPrimal::primalRay(CoinIndexedVector * rowArray)
3147	{
3148	delete [] ray_;
3149	ray_ = new double [numberColumns_];
3150	CoinZeroN(ray_, numberColumns_);
3151	int number = rowArray->getNumElements();
3152	int * index = rowArray->getIndices();
3153	double * array = rowArray->denseVector();
3154	double way = -directionIn_;
3155	int i;
3156	double zeroTolerance = 1.0e-12;
3157	if (sequenceIn_ < numberColumns_)
3158	ray_[sequenceIn_] = directionIn_;
3159	if (!rowArray->packedMode()) {
3160	for (i = 0; i < number; i++) {
3161	int iRow = index[i];
3162	int iPivot = pivotVariable_[iRow];
3163	double arrayValue = array[iRow];
3164	if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
3165	ray_[iPivot] = way * arrayValue;
3166	}
3167	} else {
3168	for (i = 0; i < number; i++) {
3169	int iRow = index[i];
3170	int iPivot = pivotVariable_[iRow];
3171	double arrayValue = array[i];
3172	if (iPivot < numberColumns_ && fabs(arrayValue) >= zeroTolerance)
3173	ray_[iPivot] = way * arrayValue;
3174	}
3175	}
3176	}
3177	/* Get next superbasic -1 if none,
3178	Normal type is 1
3179	If type is 3 then initializes sorted list
3180	if 2 uses list.
3181	*/
3182	int
3183	ClpSimplexPrimal::nextSuperBasic(int superBasicType,
3184	CoinIndexedVector * columnArray)
3185	{
3186	int returnValue = -1;
3187	bool finished = false;
3188	while (!finished) {
3189	returnValue = firstFree_;
3190	int iColumn = firstFree_ + 1;
3191	if (superBasicType > 1) {
3192	if (superBasicType > 2) {
3193	// Initialize list
3194	// Wild guess that lower bound more natural than upper
3195	int number = 0;
3196	double * work = columnArray->denseVector();
3197	int * which = columnArray->getIndices();
3198	for (iColumn = 0; iColumn < numberRows_ + numberColumns_; iColumn++) {
3199	if (!flagged(iColumn)) {
3200	if (getStatus(iColumn) == superBasic) {
3201	if (fabs(solution_[iColumn] - lower_[iColumn]) <= primalTolerance_) {
3202	solution_[iColumn] = lower_[iColumn];
3203	setStatus(iColumn, atLowerBound);
3204	} else if (fabs(solution_[iColumn] - upper_[iColumn])
3205	<= primalTolerance_) {
3206	solution_[iColumn] = upper_[iColumn];
3207	setStatus(iColumn, atUpperBound);
3208	} else if (lower_[iColumn] < -1.0e20 && upper_[iColumn] > 1.0e20) {
3209	setStatus(iColumn, isFree);
3210	break;
3211	} else if (!flagged(iColumn)) {
3212	// put ones near bounds at end after sorting
3213	work[number] = - CoinMin(0.1 * (solution_[iColumn] - lower_[iColumn]),
3214	upper_[iColumn] - solution_[iColumn]);
3215	which[number++] = iColumn;
3216	}
3217	}
3218	}
3219	}
3220	CoinSort_2(work, work + number, which);
3221	columnArray->setNumElements(number);
3222	CoinZeroN(work, number);
3223	}
3224	int * which = columnArray->getIndices();
3225	int number = columnArray->getNumElements();
3226	if (!number) {
3227	// finished
3228	iColumn = numberRows_ + numberColumns_;
3229	returnValue = -1;
3230	} else {
3231	number--;
3232	returnValue = which[number];
3233	iColumn = returnValue;
3234	columnArray->setNumElements(number);
3235	}
3236	} else {
3237	for (; iColumn < numberRows_ + numberColumns_; iColumn++) {
3238	if (!flagged(iColumn)) {
3239	if (getStatus(iColumn) == superBasic) {
3240	if (fabs(solution_[iColumn] - lower_[iColumn]) <= primalTolerance_) {
3241	solution_[iColumn] = lower_[iColumn];
3242	setStatus(iColumn, atLowerBound);
3243	} else if (fabs(solution_[iColumn] - upper_[iColumn])
3244	<= primalTolerance_) {
3245	solution_[iColumn] = upper_[iColumn];
3246	setStatus(iColumn, atUpperBound);
3247	} else if (lower_[iColumn] < -1.0e20 && upper_[iColumn] > 1.0e20) {
3248	setStatus(iColumn, isFree);
3249	break;
3250	} else {
3251	break;
3252	}
3253	}
3254	}
3255	}
3256	}
3257	firstFree_ = iColumn;
3258	finished = true;
3259	if (firstFree_ == numberRows_ + numberColumns_)
3260	firstFree_ = -1;
3261	if (returnValue >= 0 && getStatus(returnValue) != superBasic && getStatus(returnValue) != isFree)
3262	finished = false; // somehow picked up odd one
3263	}
3264	return returnValue;
3265	}
3266	void
3267	ClpSimplexPrimal::clearAll()
3268	{
3269	// Clean up any gub stuff
3270	matrix_->extendUpdated(this, rowArray_[1], 1);
3271	int number = rowArray_[1]->getNumElements();
3272	int * which = rowArray_[1]->getIndices();
3273
3274	int iIndex;
3275	for (iIndex = 0; iIndex < number; iIndex++) {
3276
3277	int iRow = which[iIndex];
3278	clearActive(iRow);
3279	}
3280	rowArray_[1]->clear();
3281	// make sure any gub sets are clean
3282	matrix_->generalExpanded(this, 11, sequenceIn_);
3283
3284	}
3285	// Sort of lexicographic resolve
3286	int
3287	ClpSimplexPrimal::lexSolve()
3288	{
3289	algorithm_ = +1;
3290	//specialOptions_ \|= 4;
3291
3292	// save data
3293	ClpDataSave data = saveData();
3294	matrix_->refresh(this); // make sure matrix okay
3295
3296	// Save so can see if doing after dual
3297	int initialStatus = problemStatus_;
3298	int initialIterations = numberIterations_;
3299	int initialNegDjs = -1;
3300	// initialize - maybe values pass and algorithm_ is +1
3301	int ifValuesPass = 0;
3302	#if 0
3303	// if so - put in any superbasic costed slacks
3304	// Start can skip some things in transposeTimes
3305	specialOptions_ \|= 131072;
3306	if (ifValuesPass && specialOptions_ < 0x01000000) {
3307	// Get column copy
3308	const CoinPackedMatrix * columnCopy = matrix();
3309	const int * row = columnCopy->getIndices();
3310	const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
3311	const int * columnLength = columnCopy->getVectorLengths();
3312	//const double * element = columnCopy->getElements();
3313	int n = 0;
3314	for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
3315	if (columnLength[iColumn] == 1) {
3316	Status status = getColumnStatus(iColumn);
3317	if (status != basic && status != isFree) {
3318	double value = columnActivity_[iColumn];
3319	if (fabs(value - columnLower_[iColumn]) > primalTolerance_ &&
3320	fabs(value - columnUpper_[iColumn]) > primalTolerance_) {
3321	int iRow = row[columnStart[iColumn]];
3322	if (getRowStatus(iRow) == basic) {
3323	setRowStatus(iRow, superBasic);
3324	setColumnStatus(iColumn, basic);
3325	n++;
3326	}
3327	}
3328	}
3329	}
3330	}
3331	printf("%d costed slacks put in basis\n", n);
3332	}
3333	#endif
3334	double * originalCost = NULL;
3335	double * originalLower = NULL;
3336	double * originalUpper = NULL;
3337	if (!startup(0, 0)) {
3338
3339	// Set average theta
3340	nonLinearCost_->setAverageTheta(1.0e3);
3341	int lastCleaned = 0; // last time objective or bounds cleaned up
3342
3343	// Say no pivot has occurred (for steepest edge and updates)
3344	pivotRow_ = -2;
3345
3346	// This says whether to restore things etc
3347	int factorType = 0;
3348	if (problemStatus_ < 0 && perturbation_ < 100) {
3349	perturb(0);
3350	// Can't get here if values pass
3351	assert (!ifValuesPass);
3352	gutsOfSolution(NULL, NULL);
3353	if (handler_->logLevel() > 2) {
3354	handler_->message(CLP_SIMPLEX_STATUS, messages_)
3355	<< numberIterations_ << objectiveValue();
3356	handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3357	<< sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3358	handler_->printing(sumDualInfeasibilities_ > 0.0)
3359	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
3360	handler_->printing(numberDualInfeasibilitiesWithoutFree_
3361	< numberDualInfeasibilities_)
3362	<< numberDualInfeasibilitiesWithoutFree_;
3363	handler_->message() << CoinMessageEol;
3364	}
3365	}
3366	ClpSimplex * saveModel = NULL;
3367	int stopSprint = -1;
3368	int sprintPass = 0;
3369	int reasonableSprintIteration = 0;
3370	int lastSprintIteration = 0;
3371	double lastObjectiveValue = COIN_DBL_MAX;
3372	// Start check for cycles
3373	progress_.fillFromModel(this);
3374	progress_.startCheck();
3375	/*
3376	Status of problem:
3377	0 - optimal
3378	1 - infeasible
3379	2 - unbounded
3380	-1 - iterating
3381	-2 - factorization wanted
3382	-3 - redo checking without factorization
3383	-4 - looks infeasible
3384	-5 - looks unbounded
3385	*/
3386	originalCost = CoinCopyOfArray(cost_, numberColumns_ + numberRows_);
3387	originalLower = CoinCopyOfArray(lower_, numberColumns_ + numberRows_);
3388	originalUpper = CoinCopyOfArray(upper_, numberColumns_ + numberRows_);
3389	while (problemStatus_ < 0) {
3390	int iRow, iColumn;
3391	// clear
3392	for (iRow = 0; iRow < 4; iRow++) {
3393	rowArray_[iRow]->clear();
3394	}
3395
3396	for (iColumn = 0; iColumn < 2; iColumn++) {
3397	columnArray_[iColumn]->clear();
3398	}
3399
3400	// give matrix (and model costs and bounds a chance to be
3401	// refreshed (normally null)
3402	matrix_->refresh(this);
3403	// If getting nowhere - why not give it a kick
3404	#if 1
3405	if (perturbation_ < 101 && numberIterations_ > 2 * (numberRows_ + numberColumns_) && (specialOptions_ & 4) == 0
3406	&& initialStatus != 10) {
3407	perturb(1);
3408	matrix_->rhsOffset(this, true, false);
3409	}
3410	#endif
3411	// If we have done no iterations - special
3412	if (lastGoodIteration_ == numberIterations_ && factorType)
3413	factorType = 3;
3414	if (saveModel) {
3415	// Doing sprint
3416	if (sequenceIn_ < 0 \|\| numberIterations_ >= stopSprint) {
3417	problemStatus_ = -1;
3418	originalModel(saveModel);
3419	saveModel = NULL;
3420	if (sequenceIn_ < 0 && numberIterations_ < reasonableSprintIteration &&
3421	sprintPass > 100)
3422	primalColumnPivot_->switchOffSprint();
3423	//lastSprintIteration=numberIterations_;
3424	printf("End small model\n");
3425	}
3426	}
3427
3428	// may factorize, checks if problem finished
3429	statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3430	if (initialStatus == 10) {
3431	// cleanup phase
3432	if(initialIterations != numberIterations_) {
3433	if (numberDualInfeasibilities_ > 10000 && numberDualInfeasibilities_ > 10 * initialNegDjs) {
3434	// getting worse - try perturbing
3435	if (perturbation_ < 101 && (specialOptions_ & 4) == 0) {
3436	perturb(1);
3437	matrix_->rhsOffset(this, true, false);
3438	statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3439	}
3440	}
3441	} else {
3442	// save number of negative djs
3443	if (!numberPrimalInfeasibilities_)
3444	initialNegDjs = numberDualInfeasibilities_;
3445	// make sure weight won't be changed
3446	if (infeasibilityCost_ == 1.0e10)
3447	infeasibilityCost_ = 1.000001e10;
3448	}
3449	}
3450	// See if sprint says redo because of problems
3451	if (numberDualInfeasibilities_ == -776) {
3452	// Need new set of variables
3453	problemStatus_ = -1;
3454	originalModel(saveModel);
3455	saveModel = NULL;
3456	//lastSprintIteration=numberIterations_;
3457	printf("End small model after\n");
3458	statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3459	}
3460	int numberSprintIterations = 0;
3461	int numberSprintColumns = primalColumnPivot_->numberSprintColumns(numberSprintIterations);
3462	if (problemStatus_ == 777) {
3463	// problems so do one pass with normal
3464	problemStatus_ = -1;
3465	originalModel(saveModel);
3466	saveModel = NULL;
3467	// Skip factorization
3468	//statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
3469	statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true, ifValuesPass, saveModel);
3470	} else if (problemStatus_ < 0 && !saveModel && numberSprintColumns && firstFree_ < 0) {
3471	int numberSort = 0;
3472	int numberFixed = 0;
3473	int numberBasic = 0;
3474	reasonableSprintIteration = numberIterations_ + 100;
3475	int * whichColumns = new int[numberColumns_];
3476	double * weight = new double[numberColumns_];
3477	int numberNegative = 0;
3478	double sumNegative = 0.0;
3479	// now massage weight so all basic in plus good djs
3480	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3481	double dj = dj_[iColumn];
3482	switch(getColumnStatus(iColumn)) {
3483
3484	case basic:
3485	dj = -1.0e50;
3486	numberBasic++;
3487	break;
3488	case atUpperBound:
3489	dj = -dj;
3490	break;
3491	case isFixed:
3492	dj = 1.0e50;
3493	numberFixed++;
3494	break;
3495	case atLowerBound:
3496	dj = dj;
3497	break;
3498	case isFree:
3499	dj = -100.0 * fabs(dj);
3500	break;
3501	case superBasic:
3502	dj = -100.0 * fabs(dj);
3503	break;
3504	}
3505	if (dj < -dualTolerance_ && dj > -1.0e50) {
3506	numberNegative++;
3507	sumNegative -= dj;
3508	}
3509	weight[iColumn] = dj;
3510	whichColumns[iColumn] = iColumn;
3511	}
3512	handler_->message(CLP_SPRINT, messages_)
3513	<< sprintPass << numberIterations_ - lastSprintIteration << objectiveValue() << sumNegative
3514	<< numberNegative
3515	<< CoinMessageEol;
3516	sprintPass++;
3517	lastSprintIteration = numberIterations_;
3518	if (objectiveValue()*optimizationDirection_ > lastObjectiveValue - 1.0e-7 && sprintPass > 5) {
3519	// switch off
3520	printf("Switching off sprint\n");
3521	primalColumnPivot_->switchOffSprint();
3522	} else {
3523	lastObjectiveValue = objectiveValue() * optimizationDirection_;
3524	// sort
3525	CoinSort_2(weight, weight + numberColumns_, whichColumns);
3526	numberSort = CoinMin(numberColumns_ - numberFixed, numberBasic + numberSprintColumns);
3527	// Sort to make consistent ?
3528	std::sort(whichColumns, whichColumns + numberSort);
3529	saveModel = new ClpSimplex(this, numberSort, whichColumns);
3530	delete [] whichColumns;
3531	delete [] weight;
3532	// Skip factorization
3533	//statusOfProblemInPrimal(lastCleaned,factorType,&progress_,false,saveModel);
3534	//statusOfProblemInPrimal(lastCleaned,factorType,&progress_,true,saveModel);
3535	stopSprint = numberIterations_ + numberSprintIterations;
3536	printf("Sprint with %d columns for %d iterations\n",
3537	numberSprintColumns, numberSprintIterations);
3538	}
3539	}
3540
3541	// Say good factorization
3542	factorType = 1;
3543
3544	// Say no pivot has occurred (for steepest edge and updates)
3545	pivotRow_ = -2;
3546
3547	// exit if victory declared
3548	if (problemStatus_ >= 0) {
3549	if (originalCost) {
3550	// find number nonbasic with zero reduced costs
3551	int numberDegen = 0;
3552	int numberTotal = numberColumns_; //+numberRows_;
3553	for (int i = 0; i < numberTotal; i++) {
3554	cost_[i] = 0.0;
3555	if (getStatus(i) == atLowerBound) {
3556	if (dj_[i] <= dualTolerance_) {
3557	cost_[i] = numberTotal - i + randomNumberGenerator_.randomDouble() * 0.5;
3558	numberDegen++;
3559	} else {
3560	// fix
3561	cost_[i] = 1.0e10; //upper_[i]=lower_[i];
3562	}
3563	} else if (getStatus(i) == atUpperBound) {
3564	if (dj_[i] >= -dualTolerance_) {
3565	cost_[i] = (numberTotal - i) + randomNumberGenerator_.randomDouble() * 0.5;
3566	numberDegen++;
3567	} else {
3568	// fix
3569	cost_[i] = -1.0e10; //lower_[i]=upper_[i];
3570	}
3571	} else if (getStatus(i) == basic) {
3572	cost_[i] = (numberTotal - i) + randomNumberGenerator_.randomDouble() * 0.5;
3573	}
3574	}
3575	problemStatus_ = -1;
3576	lastObjectiveValue = COIN_DBL_MAX;
3577	// Start check for cycles
3578	progress_.fillFromModel(this);
3579	progress_.startCheck();
3580	printf("%d degenerate after %d iterations\n", numberDegen,
3581	numberIterations_);
3582	if (!numberDegen) {
3583	CoinMemcpyN(originalCost, numberTotal, cost_);
3584	delete [] originalCost;
3585	originalCost = NULL;
3586	CoinMemcpyN(originalLower, numberTotal, lower_);
3587	delete [] originalLower;
3588	CoinMemcpyN(originalUpper, numberTotal, upper_);
3589	delete [] originalUpper;
3590	}
3591	delete nonLinearCost_;
3592	nonLinearCost_ = new ClpNonLinearCost(this);
3593	progress_.endOddState();
3594	continue;
3595	} else {
3596	printf("exiting after %d iterations\n", numberIterations_);
3597	break;
3598	}
3599	}
3600
3601	// test for maximum iterations
3602	if (hitMaximumIterations() \|\| (ifValuesPass == 2 && firstFree_ < 0)) {
3603	problemStatus_ = 3;
3604	break;
3605	}
3606
3607	if (firstFree_ < 0) {
3608	if (ifValuesPass) {
3609	// end of values pass
3610	ifValuesPass = 0;
3611	int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
3612	if (status >= 0) {
3613	problemStatus_ = 5;
3614	secondaryStatus_ = ClpEventHandler::endOfValuesPass;
3615	break;
3616	}
3617	}
3618	}
3619	// Check event
3620	{
3621	int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3622	if (status >= 0) {
3623	problemStatus_ = 5;
3624	secondaryStatus_ = ClpEventHandler::endOfFactorization;
3625	break;
3626	}
3627	}
3628	// Iterate
3629	whileIterating(ifValuesPass ? 1 : 0);
3630	}
3631	}
3632	assert (!originalCost);
3633	// if infeasible get real values
3634	//printf("XXXXY final cost %g\n",infeasibilityCost_);
3635	progress_.initialWeight_ = 0.0;
3636	if (problemStatus_ == 1 && secondaryStatus_ != 6) {
3637	infeasibilityCost_ = 0.0;
3638	createRim(1 + 4);
3639	nonLinearCost_->checkInfeasibilities(0.0);
3640	sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
3641	numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
3642	// and get good feasible duals
3643	computeDuals(NULL);
3644	}
3645	// Stop can skip some things in transposeTimes
3646	specialOptions_ &= ~131072;
3647	// clean up
3648	unflag();
3649	finish(0);
3650	restoreData(data);
3651	return problemStatus_;
3652	}

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

Download in other formats: