Context Navigation

ClpSimplex.cpp @ 2114

Last change on this file since 2114 was 2114, checked in by forrest, 5 years ago
was not totally threadsafe!
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 523.1 KB

Line
1	/* $Id: ClpSimplex.cpp 2114 2015-02-10 12:12:46Z forrest $ */
2	// Copyright (C) 2002, International Business Machines
3	// Corporation and others. All Rights Reserved.
4	// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6	//#undef NDEBUG
7
8	#include "ClpConfig.h"
9
10	#include "CoinPragma.hpp"
11	#include <math.h>
12
13	#if SLIM_CLP==2
14	#define SLIM_NOIO
15	#endif
16	#include "CoinHelperFunctions.hpp"
17	#include "CoinFloatEqual.hpp"
18	#include "ClpSimplex.hpp"
19	#include "ClpFactorization.hpp"
20	#include "ClpPackedMatrix.hpp"
21	#include "CoinIndexedVector.hpp"
22	#include "ClpDualRowDantzig.hpp"
23	#include "ClpDualRowSteepest.hpp"
24	#include "ClpPrimalColumnDantzig.hpp"
25	#include "ClpPrimalColumnSteepest.hpp"
26	#include "ClpNonLinearCost.hpp"
27	#include "ClpMessage.hpp"
28	#include "ClpEventHandler.hpp"
29	#include "ClpLinearObjective.hpp"
30	#include "ClpHelperFunctions.hpp"
31	#include "CoinModel.hpp"
32	#include "CoinLpIO.hpp"
33	#include <cfloat>
34	#if CLP_HAS_ABC
35	#include "CoinAbcCommon.hpp"
36	#endif
37
38	#include <string>
39	#include <stdio.h>
40	#include <iostream>
41	//#############################################################################
42
43	ClpSimplex::ClpSimplex (bool emptyMessages) :
44
45	ClpModel(emptyMessages),
46	bestPossibleImprovement_(0.0),
47	zeroTolerance_(1.0e-13),
48	columnPrimalSequence_(-2),
49	rowPrimalSequence_(-2),
50	bestObjectiveValue_(-COIN_DBL_MAX),
51	moreSpecialOptions_(2),
52	baseIteration_(0),
53	primalToleranceToGetOptimal_(-1.0),
54	largeValue_(1.0e15),
55	largestPrimalError_(0.0),
56	largestDualError_(0.0),
57	alphaAccuracy_(-1.0),
58	dualBound_(1.0e10),
59	alpha_(0.0),
60	theta_(0.0),
61	lowerIn_(0.0),
62	valueIn_(0.0),
63	upperIn_(-COIN_DBL_MAX),
64	dualIn_(0.0),
65	lowerOut_(-1),
66	valueOut_(-1),
67	upperOut_(-1),
68	dualOut_(-1),
69	dualTolerance_(1.0e-7),
70	primalTolerance_(1.0e-7),
71	sumDualInfeasibilities_(0.0),
72	sumPrimalInfeasibilities_(0.0),
73	infeasibilityCost_(1.0e10),
74	sumOfRelaxedDualInfeasibilities_(0.0),
75	sumOfRelaxedPrimalInfeasibilities_(0.0),
76	acceptablePivot_(1.0e-8),
77	lower_(NULL),
78	rowLowerWork_(NULL),
79	columnLowerWork_(NULL),
80	upper_(NULL),
81	rowUpperWork_(NULL),
82	columnUpperWork_(NULL),
83	cost_(NULL),
84	rowObjectiveWork_(NULL),
85	objectiveWork_(NULL),
86	sequenceIn_(-1),
87	directionIn_(-1),
88	sequenceOut_(-1),
89	directionOut_(-1),
90	pivotRow_(-1),
91	lastGoodIteration_(-100),
92	dj_(NULL),
93	rowReducedCost_(NULL),
94	reducedCostWork_(NULL),
95	solution_(NULL),
96	rowActivityWork_(NULL),
97	columnActivityWork_(NULL),
98	numberDualInfeasibilities_(0),
99	numberDualInfeasibilitiesWithoutFree_(0),
100	numberPrimalInfeasibilities_(100),
101	numberRefinements_(0),
102	pivotVariable_(NULL),
103	factorization_(NULL),
104	savedSolution_(NULL),
105	numberTimesOptimal_(0),
106	disasterArea_(NULL),
107	changeMade_(1),
108	algorithm_(0),
109	forceFactorization_(-1),
110	perturbation_(100),
111	nonLinearCost_(NULL),
112	lastBadIteration_(-999999),
113	lastFlaggedIteration_(-999999),
114	numberFake_(0),
115	numberChanged_(0),
116	progressFlag_(0),
117	firstFree_(-1),
118	numberExtraRows_(0),
119	maximumBasic_(0),
120	dontFactorizePivots_(0),
121	incomingInfeasibility_(1.0),
122	allowedInfeasibility_(10.0),
123	automaticScale_(0),
124	maximumPerturbationSize_(0),
125	perturbationArray_(NULL),
126	baseModel_(NULL)
127	#ifdef ABC_INHERIT
128	,abcSimplex_(NULL),
129	abcState_(0)
130	#endif
131	{
132	int i;
133	for (i = 0; i < 6; i++) {
134	rowArray_[i] = NULL;
135	columnArray_[i] = NULL;
136	}
137	for (i = 0; i < 4; i++) {
138	spareIntArray_[i] = 0;
139	spareDoubleArray_[i] = 0.0;
140	}
141	saveStatus_ = NULL;
142	// get an empty factorization so we can set tolerances etc
143	getEmptyFactorization();
144	// Say sparse
145	factorization_->sparseThreshold(1);
146	// say Steepest pricing
147	dualRowPivot_ = new ClpDualRowSteepest();
148	// say Steepest pricing
149	primalColumnPivot_ = new ClpPrimalColumnSteepest();
150	solveType_ = 1; // say simplex based life form
151
152	}
153
154	// Subproblem constructor
155	ClpSimplex::ClpSimplex ( const ClpModel * rhs,
156	int numberRows, const int * whichRow,
157	int numberColumns, const int * whichColumn,
158	bool dropNames, bool dropIntegers, bool fixOthers)
159	: ClpModel(rhs, numberRows, whichRow,
160	numberColumns, whichColumn, dropNames, dropIntegers),
161	bestPossibleImprovement_(0.0),
162	zeroTolerance_(1.0e-13),
163	columnPrimalSequence_(-2),
164	rowPrimalSequence_(-2),
165	bestObjectiveValue_(-COIN_DBL_MAX),
166	moreSpecialOptions_(2),
167	baseIteration_(0),
168	primalToleranceToGetOptimal_(-1.0),
169	largeValue_(1.0e15),
170	largestPrimalError_(0.0),
171	largestDualError_(0.0),
172	alphaAccuracy_(-1.0),
173	dualBound_(1.0e10),
174	alpha_(0.0),
175	theta_(0.0),
176	lowerIn_(0.0),
177	valueIn_(0.0),
178	upperIn_(-COIN_DBL_MAX),
179	dualIn_(0.0),
180	lowerOut_(-1),
181	valueOut_(-1),
182	upperOut_(-1),
183	dualOut_(-1),
184	dualTolerance_(1.0e-7),
185	primalTolerance_(1.0e-7),
186	sumDualInfeasibilities_(0.0),
187	sumPrimalInfeasibilities_(0.0),
188	infeasibilityCost_(1.0e10),
189	sumOfRelaxedDualInfeasibilities_(0.0),
190	sumOfRelaxedPrimalInfeasibilities_(0.0),
191	acceptablePivot_(1.0e-8),
192	lower_(NULL),
193	rowLowerWork_(NULL),
194	columnLowerWork_(NULL),
195	upper_(NULL),
196	rowUpperWork_(NULL),
197	columnUpperWork_(NULL),
198	cost_(NULL),
199	rowObjectiveWork_(NULL),
200	objectiveWork_(NULL),
201	sequenceIn_(-1),
202	directionIn_(-1),
203	sequenceOut_(-1),
204	directionOut_(-1),
205	pivotRow_(-1),
206	lastGoodIteration_(-100),
207	dj_(NULL),
208	rowReducedCost_(NULL),
209	reducedCostWork_(NULL),
210	solution_(NULL),
211	rowActivityWork_(NULL),
212	columnActivityWork_(NULL),
213	numberDualInfeasibilities_(0),
214	numberDualInfeasibilitiesWithoutFree_(0),
215	numberPrimalInfeasibilities_(100),
216	numberRefinements_(0),
217	pivotVariable_(NULL),
218	factorization_(NULL),
219	savedSolution_(NULL),
220	numberTimesOptimal_(0),
221	disasterArea_(NULL),
222	changeMade_(1),
223	algorithm_(0),
224	forceFactorization_(-1),
225	perturbation_(100),
226	nonLinearCost_(NULL),
227	lastBadIteration_(-999999),
228	lastFlaggedIteration_(-999999),
229	numberFake_(0),
230	numberChanged_(0),
231	progressFlag_(0),
232	firstFree_(-1),
233	numberExtraRows_(0),
234	maximumBasic_(0),
235	dontFactorizePivots_(0),
236	incomingInfeasibility_(1.0),
237	allowedInfeasibility_(10.0),
238	automaticScale_(0),
239	maximumPerturbationSize_(0),
240	perturbationArray_(NULL),
241	baseModel_(NULL)
242	#ifdef ABC_INHERIT
243	,abcSimplex_(NULL),
244	abcState_(0)
245	#endif
246	{
247	int i;
248	for (i = 0; i < 6; i++) {
249	rowArray_[i] = NULL;
250	columnArray_[i] = NULL;
251	}
252	for (i = 0; i < 4; i++) {
253	spareIntArray_[i] = 0;
254	spareDoubleArray_[i] = 0.0;
255	}
256	saveStatus_ = NULL;
257	// get an empty factorization so we can set tolerances etc
258	getEmptyFactorization();
259	// say Steepest pricing
260	dualRowPivot_ = new ClpDualRowSteepest();
261	// say Steepest pricing
262	primalColumnPivot_ = new ClpPrimalColumnSteepest();
263	solveType_ = 1; // say simplex based life form
264	if (fixOthers) {
265	int numberOtherColumns = rhs->numberColumns();
266	int numberOtherRows = rhs->numberRows();
267	double * solution = new double [numberOtherColumns];
268	CoinZeroN(solution, numberOtherColumns);
269	int i;
270	for (i = 0; i < numberColumns; i++) {
271	int iColumn = whichColumn[i];
272	if (solution[iColumn])
273	fixOthers = false; // duplicates
274	solution[iColumn] = 1.0;
275	}
276	if (fixOthers) {
277	const double * otherSolution = rhs->primalColumnSolution();
278	const double * objective = rhs->objective();
279	double offset = 0.0;
280	for (i = 0; i < numberOtherColumns; i++) {
281	if (solution[i]) {
282	solution[i] = 0.0; // in
283	} else {
284	solution[i] = otherSolution[i];
285	offset += objective[i] * otherSolution[i];
286	}
287	}
288	double * rhsModification = new double [numberOtherRows];
289	CoinZeroN(rhsModification, numberOtherRows);
290	rhs->matrix()->times(solution, rhsModification) ;
291	for ( i = 0; i < numberRows; i++) {
292	int iRow = whichRow[i];
293	if (rowLower_[i] > -1.0e20)
294	rowLower_[i] -= rhsModification[iRow];
295	if (rowUpper_[i] < 1.0e20)
296	rowUpper_[i] -= rhsModification[iRow];
297	}
298	delete [] rhsModification;
299	setObjectiveOffset(rhs->objectiveOffset() - offset);
300	// And set objective value to match
301	setObjectiveValue(rhs->objectiveValue());
302	}
303	delete [] solution;
304	}
305	}
306	// Subproblem constructor
307	ClpSimplex::ClpSimplex ( const ClpSimplex * rhs,
308	int numberRows, const int * whichRow,
309	int numberColumns, const int * whichColumn,
310	bool dropNames, bool dropIntegers, bool fixOthers)
311	: ClpModel(rhs, numberRows, whichRow,
312	numberColumns, whichColumn, dropNames, dropIntegers),
313	bestPossibleImprovement_(0.0),
314	zeroTolerance_(1.0e-13),
315	columnPrimalSequence_(-2),
316	rowPrimalSequence_(-2),
317	bestObjectiveValue_(-COIN_DBL_MAX),
318	moreSpecialOptions_(2),
319	baseIteration_(0),
320	primalToleranceToGetOptimal_(-1.0),
321	largeValue_(1.0e15),
322	largestPrimalError_(0.0),
323	largestDualError_(0.0),
324	alphaAccuracy_(-1.0),
325	dualBound_(1.0e10),
326	alpha_(0.0),
327	theta_(0.0),
328	lowerIn_(0.0),
329	valueIn_(0.0),
330	upperIn_(-COIN_DBL_MAX),
331	dualIn_(0.0),
332	lowerOut_(-1),
333	valueOut_(-1),
334	upperOut_(-1),
335	dualOut_(-1),
336	dualTolerance_(rhs->dualTolerance_),
337	primalTolerance_(rhs->primalTolerance_),
338	sumDualInfeasibilities_(0.0),
339	sumPrimalInfeasibilities_(0.0),
340	infeasibilityCost_(1.0e10),
341	sumOfRelaxedDualInfeasibilities_(0.0),
342	sumOfRelaxedPrimalInfeasibilities_(0.0),
343	acceptablePivot_(1.0e-8),
344	lower_(NULL),
345	rowLowerWork_(NULL),
346	columnLowerWork_(NULL),
347	upper_(NULL),
348	rowUpperWork_(NULL),
349	columnUpperWork_(NULL),
350	cost_(NULL),
351	rowObjectiveWork_(NULL),
352	objectiveWork_(NULL),
353	sequenceIn_(-1),
354	directionIn_(-1),
355	sequenceOut_(-1),
356	directionOut_(-1),
357	pivotRow_(-1),
358	lastGoodIteration_(-100),
359	dj_(NULL),
360	rowReducedCost_(NULL),
361	reducedCostWork_(NULL),
362	solution_(NULL),
363	rowActivityWork_(NULL),
364	columnActivityWork_(NULL),
365	numberDualInfeasibilities_(0),
366	numberDualInfeasibilitiesWithoutFree_(0),
367	numberPrimalInfeasibilities_(100),
368	numberRefinements_(0),
369	pivotVariable_(NULL),
370	factorization_(NULL),
371	savedSolution_(NULL),
372	numberTimesOptimal_(0),
373	disasterArea_(NULL),
374	changeMade_(1),
375	algorithm_(0),
376	forceFactorization_(-1),
377	perturbation_(100),
378	nonLinearCost_(NULL),
379	lastBadIteration_(-999999),
380	lastFlaggedIteration_(-999999),
381	numberFake_(0),
382	numberChanged_(0),
383	progressFlag_(0),
384	firstFree_(-1),
385	numberExtraRows_(0),
386	maximumBasic_(0),
387	dontFactorizePivots_(0),
388	incomingInfeasibility_(1.0),
389	allowedInfeasibility_(10.0),
390	automaticScale_(0),
391	maximumPerturbationSize_(0),
392	perturbationArray_(NULL),
393	baseModel_(NULL)
394	#ifdef ABC_INHERIT
395	,abcSimplex_(NULL),
396	abcState_(rhs->abcState_)
397	#endif
398	{
399	int i;
400	for (i = 0; i < 6; i++) {
401	rowArray_[i] = NULL;
402	columnArray_[i] = NULL;
403	}
404	for (i = 0; i < 4; i++) {
405	spareIntArray_[i] = 0;
406	spareDoubleArray_[i] = 0.0;
407	}
408	saveStatus_ = NULL;
409	factorization_ = new ClpFactorization(*rhs->factorization_, -numberRows_);
410	//factorization_ = new ClpFactorization(*rhs->factorization_,
411	// rhs->factorization_->goDenseThreshold());
412	ClpDualRowDantzig * pivot =
413	dynamic_cast< ClpDualRowDantzig*>(rhs->dualRowPivot_);
414	// say Steepest pricing
415	if (!pivot)
416	dualRowPivot_ = new ClpDualRowSteepest();
417	else
418	dualRowPivot_ = new ClpDualRowDantzig();
419	// say Steepest pricing
420	primalColumnPivot_ = new ClpPrimalColumnSteepest();
421	solveType_ = 1; // say simplex based life form
422	if (fixOthers) {
423	int numberOtherColumns = rhs->numberColumns();
424	int numberOtherRows = rhs->numberRows();
425	double * solution = new double [numberOtherColumns];
426	CoinZeroN(solution, numberOtherColumns);
427	int i;
428	for (i = 0; i < numberColumns; i++) {
429	int iColumn = whichColumn[i];
430	if (solution[iColumn])
431	fixOthers = false; // duplicates
432	solution[iColumn] = 1.0;
433	}
434	if (fixOthers) {
435	const double * otherSolution = rhs->primalColumnSolution();
436	const double * objective = rhs->objective();
437	double offset = 0.0;
438	for (i = 0; i < numberOtherColumns; i++) {
439	if (solution[i]) {
440	solution[i] = 0.0; // in
441	} else {
442	solution[i] = otherSolution[i];
443	offset += objective[i] * otherSolution[i];
444	}
445	}
446	double * rhsModification = new double [numberOtherRows];
447	CoinZeroN(rhsModification, numberOtherRows);
448	rhs->matrix()->times(solution, rhsModification) ;
449	for ( i = 0; i < numberRows; i++) {
450	int iRow = whichRow[i];
451	if (rowLower_[i] > -1.0e20)
452	rowLower_[i] -= rhsModification[iRow];
453	if (rowUpper_[i] < 1.0e20)
454	rowUpper_[i] -= rhsModification[iRow];
455	}
456	delete [] rhsModification;
457	setObjectiveOffset(rhs->objectiveOffset() - offset);
458	// And set objective value to match
459	setObjectiveValue(rhs->objectiveValue());
460	}
461	delete [] solution;
462	}
463	if (rhs->maximumPerturbationSize_) {
464	maximumPerturbationSize_ = 2 * numberColumns;
465	perturbationArray_ = new double [maximumPerturbationSize_];
466	for (i = 0; i < numberColumns; i++) {
467	int iColumn = whichColumn[i];
468	perturbationArray_[2i] = rhs->perturbationArray_[2iColumn];
469	perturbationArray_[2i+1] = rhs->perturbationArray_[2iColumn+1];
470	}
471	}
472	}
473	// Puts solution back small model
474	void
475	ClpSimplex::getbackSolution(const ClpSimplex & smallModel, const int * whichRow, const int * whichColumn)
476	{
477	setSumDualInfeasibilities(smallModel.sumDualInfeasibilities());
478	setNumberDualInfeasibilities(smallModel.numberDualInfeasibilities());
479	setSumPrimalInfeasibilities(smallModel.sumPrimalInfeasibilities());
480	setNumberPrimalInfeasibilities(smallModel.numberPrimalInfeasibilities());
481	setNumberIterations(smallModel.numberIterations());
482	setProblemStatus(smallModel.status());
483	setObjectiveValue(smallModel.objectiveValue());
484	const double * solution2 = smallModel.primalColumnSolution();
485	int i;
486	int numberRows2 = smallModel.numberRows();
487	int numberColumns2 = smallModel.numberColumns();
488	const double * dj2 = smallModel.dualColumnSolution();
489	for ( i = 0; i < numberColumns2; i++) {
490	int iColumn = whichColumn[i];
491	columnActivity_[iColumn] = solution2[i];
492	reducedCost_[iColumn] = dj2[i];
493	setStatus(iColumn, smallModel.getStatus(i));
494	}
495	const double * dual2 = smallModel.dualRowSolution();
496	memset(dual_, 0, numberRows_ * sizeof(double));
497	for (i = 0; i < numberRows2; i++) {
498	int iRow = whichRow[i];
499	setRowStatus(iRow, smallModel.getRowStatus(i));
500	dual_[iRow] = dual2[i];
501	}
502	CoinZeroN(rowActivity_, numberRows_);
503	#if 0
504	if (!problemStatus_) {
505	ClpDisjointCopyN(smallModel.objective(), smallModel.numberColumns_, smallModel.reducedCost_);
506	smallModel.matrix_->transposeTimes(-1.0, smallModel.dual_, smallModel.reducedCost_);
507	for (int i = 0; i < smallModel.numberColumns_; i++) {
508	if (smallModel.getColumnStatus(i) == basic)
509	assert (fabs(smallModel.reducedCost_[i]) < 1.0e-5);
510	}
511	ClpDisjointCopyN(objective(), numberColumns_, reducedCost_);
512	matrix_->transposeTimes(-1.0, dual_, reducedCost_);
513	for (int i = 0; i < numberColumns_; i++) {
514	if (getColumnStatus(i) == basic)
515	assert (fabs(reducedCost_[i]) < 1.0e-5);
516	}
517	}
518	#endif
519	matrix()->times(columnActivity_, rowActivity_) ;
520	}
521
522	//-----------------------------------------------------------------------------
523
524	ClpSimplex::~ClpSimplex ()
525	{
526	setPersistenceFlag(0);
527	gutsOfDelete(0);
528	delete nonLinearCost_;
529	}
530	//#############################################################################
531	void ClpSimplex::setLargeValue( double value)
532	{
533	if (value > 0.0 && value < COIN_DBL_MAX)
534	largeValue_ = value;
535	}
536	int
537	ClpSimplex::gutsOfSolution ( double * givenDuals,
538	const double * givenPrimals,
539	bool valuesPass)
540	{
541
542
543	// if values pass, save values of basic variables
544	double * save = NULL;
545	double oldValue = 0.0;
546	double oldLargestPrimalError=largestPrimalError_;
547	double oldLargestDualError=largestDualError_;
548	if (valuesPass) {
549	assert(algorithm_ > 0); // only primal at present
550	assert(nonLinearCost_);
551	int iRow;
552	checkPrimalSolution( rowActivityWork_, columnActivityWork_);
553	// get correct bounds on all variables
554	nonLinearCost_->checkInfeasibilities(primalTolerance_);
555	oldValue = nonLinearCost_->largestInfeasibility();
556	save = new double[numberRows_];
557	for (iRow = 0; iRow < numberRows_; iRow++) {
558	int iPivot = pivotVariable_[iRow];
559	save[iRow] = solution_[iPivot];
560	}
561	}
562	// do work
563	computePrimals(rowActivityWork_, columnActivityWork_);
564	// If necessary - override results
565	if (givenPrimals) {
566	CoinMemcpyN(givenPrimals, numberColumns_, columnActivityWork_);
567	memset(rowActivityWork_, 0, numberRows_ * sizeof(double));
568	times(-1.0, columnActivityWork_, rowActivityWork_);
569	}
570	double objectiveModification = 0.0;
571	if (algorithm_ > 0 && nonLinearCost_ != NULL) {
572	// primal algorithm
573	// get correct bounds on all variables
574	// If 4 bit set - Force outgoing variables to exact bound (primal)
575	if ((specialOptions_ & 4) == 0)
576	nonLinearCost_->checkInfeasibilities(primalTolerance_);
577	else
578	nonLinearCost_->checkInfeasibilities(0.0);
579	objectiveModification += nonLinearCost_->changeInCost();
580	if (nonLinearCost_->numberInfeasibilities())
581	if (handler_->detail(CLP_SIMPLEX_NONLINEAR, messages_) < 100) {
582	handler_->message(CLP_SIMPLEX_NONLINEAR, messages_)
583	<< nonLinearCost_->changeInCost()
584	<< nonLinearCost_->numberInfeasibilities()
585	<< CoinMessageEol;
586	}
587	}
588	if (valuesPass) {
589	double badInfeasibility = nonLinearCost_->largestInfeasibility();
590	#ifdef CLP_DEBUG
591	std::cout << "Largest given infeasibility " << oldValue
592	<< " now " << nonLinearCost_->largestInfeasibility() << std::endl;
593	#endif
594	int numberOut = 0;
595	// But may be very large rhs etc
596	double useError = CoinMin(largestPrimalError_,
597	1.0e5 / maximumAbsElement(solution_, numberRows_ + numberColumns_));
598	if ((oldValue < incomingInfeasibility_ \|\| badInfeasibility >
599	(CoinMax(10.0 * allowedInfeasibility_, 100.0 * oldValue)))
600	&& (badInfeasibility > CoinMax(incomingInfeasibility_, allowedInfeasibility_) \|\|
601	useError > 1.0e-3)) {
602	if (algorithm_>1) {
603	// nonlinear
604	//printf("Original largest infeas %g, now %g, primalError %g\n",
605	// oldValue,nonLinearCost_->largestInfeasibility(),
606	// largestPrimalError_);
607	//printf("going to all slack\n");
608	allSlackBasis(true);
609	CoinIotaN(pivotVariable_, numberRows_, numberColumns_);
610	return 1;
611	}
612	//printf("Original largest infeas %g, now %g, primalError %g\n",
613	// oldValue,nonLinearCost_->largestInfeasibility(),
614	// largestPrimalError_);
615	// throw out up to 1000 structurals
616	int maxOut = (allowedInfeasibility_==10.0) ? 1000 : 100;
617	int iRow;
618	int * sort = new int[numberRows_];
619	// first put back solution and store difference
620	for (iRow = 0; iRow < numberRows_; iRow++) {
621	int iPivot = pivotVariable_[iRow];
622	double difference = fabs(solution_[iPivot] - save[iRow]);
623	solution_[iPivot] = save[iRow];
624	save[iRow] = difference;
625	}
626	int numberBasic = 0;
627	for (iRow = 0; iRow < numberRows_; iRow++) {
628	int iPivot = pivotVariable_[iRow];
629
630	if (iPivot < numberColumns_) {
631	// column
632	double difference = save[iRow];
633	if (difference > 1.0e-4) {
634	sort[numberOut] = iRow;
635	save[numberOut++] = -difference;
636	if (getStatus(iPivot) == basic)
637	numberBasic++;
638	}
639	}
640	}
641	if (!numberBasic) {
642	//printf("no errors on basic - going to all slack - numberOut %d\n",numberOut);
643	#if 0
644	allSlackBasis(true);
645	CoinIotaN(pivotVariable_, numberRows_, numberColumns_);
646	#else
647	// allow
648	numberOut = 0;
649	#endif
650	}
651	CoinSort_2(save, save + numberOut, sort);
652	numberOut = CoinMin(maxOut, numberOut);
653	for (iRow = 0; iRow < numberOut; iRow++) {
654	int jRow = sort[iRow];
655	int iColumn = pivotVariable_[jRow];
656	setColumnStatus(iColumn, superBasic);
657	setRowStatus(jRow, basic);
658	pivotVariable_[jRow] = jRow + numberColumns_;
659	if (fabs(solution_[iColumn]) > 1.0e10) {
660	if (upper_[iColumn] < 0.0) {
661	solution_[iColumn] = upper_[iColumn];
662	} else if (lower_[iColumn] > 0.0) {
663	solution_[iColumn] = lower_[iColumn];
664	} else {
665	solution_[iColumn] = 0.0;
666	}
667	}
668	}
669	delete [] sort;
670	}
671	delete [] save;
672	if (numberOut)
673	return numberOut;
674	}
675	if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
676	//printf("trying feas pump\n");
677	const char * integerType = integerInformation();
678	assert (integerType);
679	assert (perturbationArray_);
680	CoinZeroN(cost_, numberRows_ + numberColumns_);
681	for (int i = 0; i < numberRows_ - numberRows_; i++) {
682	int iSequence = pivotVariable_[i];
683	if (iSequence < numberColumns_ && integerType[iSequence]) {
684	double lower = lower_[iSequence];
685	double upper = upper_[iSequence];
686	double value = solution_[iSequence];
687	if (value >= lower - primalTolerance_ &&
688	value <= upper + primalTolerance_) {
689	double sign;
690	if (value - lower < upper - value)
691	sign = 1.0;
692	else
693	sign = -1.0;
694	cost_[iSequence] = sign * perturbationArray_[iSequence];
695	}
696	}
697	}
698	}
699	#if CAN_HAVE_ZERO_OBJ>1
700	if ((specialOptions_&2097152)==0) {
701	#endif
702	computeDuals(givenDuals);
703	if ((moreSpecialOptions_ & 128) != 0 && !numberIterations_) {
704	const char * integerType = integerInformation();
705	// Need to do columns and rows to stay dual feasible
706	for (int iSequence = 0; iSequence < numberColumns_; iSequence++) {
707	if (integerType[iSequence] && getStatus(iSequence) != basic) {
708	double djValue = dj_[iSequence];
709	double change = 0.0;
710	if (getStatus(iSequence) == atLowerBound)
711	change = CoinMax(-djValue, 10.0 * perturbationArray_[iSequence]);
712	else if (getStatus(iSequence) == atUpperBound)
713	change = CoinMin(-djValue, -10.0 * perturbationArray_[iSequence]);
714	cost_[iSequence] = change;
715	dj_[iSequence] += change;
716	}
717	}
718	}
719
720	// now check solutions
721	//checkPrimalSolution( rowActivityWork_, columnActivityWork_);
722	//checkDualSolution();
723	checkBothSolutions();
724	objectiveValue_ += objectiveModification / (objectiveScale_ * rhsScale_);
725	#if CAN_HAVE_ZERO_OBJ>1
726	} else {
727	checkPrimalSolution( rowActivityWork_, columnActivityWork_);
728	#ifndef COIN_REUSE_RANDOM
729	memset(dj_,0,(numberRows_+numberColumns_)*sizeof(double));
730	#else
731	for (int iSequence=0;iSequence<numberRows_+numberColumns_;iSequence++) {
732	double value;
733	switch (getStatus(iSequence)) {
734	case atLowerBound:
735	value=1.0e-9*(1.0+CoinDrand48());
736	break;
737	case atUpperBound:
738	value=-1.0e-9*(1.0+CoinDrand48());
739	break;
740	default:
741	value=0.0;
742	break;
743	}
744	}
745	#endif
746	objectiveValue_=0.0;
747	}
748	#endif
749	if (handler_->logLevel() > 3 \|\| (largestPrimalError_ > 1.0e-2 \|\|
750	largestDualError_ > 1.0e-2))
751	handler_->message(CLP_SIMPLEX_ACCURACY, messages_)
752	<< largestPrimalError_
753	<< largestDualError_
754	<< CoinMessageEol;
755	if (largestPrimalError_ > 1.0e-1 && numberRows_ > 100 && numberIterations_) {
756	// Change factorization tolerance
757	if (factorization_->zeroTolerance() > 1.0e-18)
758	factorization_->zeroTolerance(1.0e-18);
759	}
760	int returnCode=0;
761	bool notChanged=true;
762	if (numberIterations_ && (forceFactorization_ > 2 \|\| forceFactorization_<0 \|\| factorization_->pivotTolerance()<0.9899999999) &&
763	(oldLargestDualError\|\|oldLargestPrimalError)) {
764	double useOldDualError=oldLargestDualError;
765	double useDualError=largestDualError_;
766	if (algorithm_>0&&nonLinearCost_&&
767	nonLinearCost_->sumInfeasibilities()) {
768	double factor=CoinMax(1.0,CoinMin(1.0e3,infeasibilityCost_*1.0e-6));
769	useOldDualError /= factor;
770	useDualError /= factor;
771	}
772	if ((largestPrimalError_>1.0e3&&
773	oldLargestPrimalError*1.0e2<largestPrimalError_)\|\|
774	(useDualError>1.0e3&&
775	useOldDualError*1.0e2<useDualError)) {
776	double pivotTolerance = factorization_->pivotTolerance();
777	double factor=(largestPrimalError_>1.0e10\|\|largestDualError_>1.0e10)
778	? 2.0 : 1.2;
779	if (pivotTolerance<0.1)
780	factorization_->pivotTolerance(0.1);
781	else if (pivotTolerance<0.98999999)
782	factorization_->pivotTolerance(CoinMin(0.99,pivotTolerance*factor));
783	notChanged=pivotTolerance==factorization_->pivotTolerance();
784	#ifdef CLP_USEFUL_PRINTOUT
785	if (pivotTolerance<0.9899999) {
786	printf("Changing pivot tolerance from %g to %g and backtracking\n",
787	pivotTolerance,factorization_->pivotTolerance());
788	}
789	printf("because old,new primal error %g,%g - dual %g,%g pivot_tol %g\n",
790	oldLargestPrimalError,largestPrimalError_,
791	oldLargestDualError,largestDualError_,
792	pivotTolerance);
793	#endif
794	if (pivotTolerance<0.9899999) {
795	largestPrimalError_=0.0;
796	largestDualError_=0.0;
797	returnCode=1;
798	}
799	}
800	}
801	if (progress_.iterationNumber_[0]>0&&
802	progress_.iterationNumber_[CLP_PROGRESS-1]
803	-progress_.iterationNumber_[0]<CLP_PROGRESS*3&&
804	factorization_->pivotTolerance()<0.25&&notChanged) {
805	double pivotTolerance = factorization_->pivotTolerance();
806	factorization_->pivotTolerance(pivotTolerance*1.5);
807	#ifdef CLP_USEFUL_PRINTOUT
808	printf("Changing pivot tolerance from %g to %g - inverting too often\n",
809	pivotTolerance,factorization_->pivotTolerance());
810	#endif
811	}
812	// Switch off false values pass indicator
813	if (!valuesPass && algorithm_ > 0)
814	firstFree_ = -1;
815	if (handler_->logLevel()==63)
816	printf("end getsolution algorithm %d status %d npinf %d sum,relaxed %g,%g ndinf %d sum,relaxed %g,%g\n",
817	algorithm_,problemStatus_,
818	numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,sumOfRelaxedPrimalInfeasibilities_,
819	numberDualInfeasibilities_,sumDualInfeasibilities_,sumOfRelaxedDualInfeasibilities_);
820	if ((moreSpecialOptions_&8388608)!=0) {
821	if (algorithm_<0) {
822	bool doneFiddling=false;
823	// Optimization may make exact test iffy
824	double testTolerance=minimumPrimalTolerance_+1.0e-15;
825	while (!doneFiddling) {
826	doneFiddling=true;
827	while( !sumOfRelaxedPrimalInfeasibilities_&&
828	primalTolerance_>testTolerance) {
829	// feasible - adjust tolerance
830	double saveTolerance=primalTolerance_;
831	primalTolerance_=CoinMax(0.25*primalTolerance_,
832	minimumPrimalTolerance_);
833	printf("Resetting primal tolerance from %g to %g\n",
834	saveTolerance,primalTolerance_);
835	dblParam_[ClpPrimalTolerance]=primalTolerance_;
836	moreSpecialOptions_ &= ~8388608;
837	// redo with switch off
838	returnCode=gutsOfSolution ( givenDuals,givenPrimals,valuesPass);
839	}
840	if(primalTolerance_>testTolerance)
841	moreSpecialOptions_ \|= 8388608; // back on
842	if ((moreSpecialOptions_&8388608)!=0) {
843	assert( numberPrimalInfeasibilities_);
844	// average infeasibility
845	double average=sumPrimalInfeasibilities_/numberPrimalInfeasibilities_;
846	double minimum=COIN_DBL_MAX;
847	double averageTotal=average;
848	bool firstTime=averageInfeasibility_[0]==COIN_DBL_MAX;
849	for (int i=0;i<CLP_INFEAS_SAVE-1;i++) {
850	double value = averageInfeasibility_[i+1];
851	averageTotal+=value;
852	averageInfeasibility_[i]=value;
853	minimum=CoinMin(minimum,value);
854	}
855	averageInfeasibility_[CLP_INFEAS_SAVE-1]=average;
856	averageTotal /= CLP_INFEAS_SAVE;
857	double oldTolerance=primalTolerance_;
858	if (averageInfeasibility_[0]!=COIN_DBL_MAX) {
859	if (firstTime) {
860	primalTolerance_=CoinMin(0.1,0.1*averageTotal);
861	primalTolerance_ = CoinMin(primalTolerance_,average);
862	} else if (primalTolerance_>0.1*minimum) {
863	primalTolerance_=0.1*minimum;
864	}
865	primalTolerance_=
866	CoinMax(primalTolerance_,minimumPrimalTolerance_);
867	}
868	if (primalTolerance_!=oldTolerance) {
869	printf("Changing primal tolerance from %g to %g\n",
870	oldTolerance,primalTolerance_);
871	moreSpecialOptions_ &= ~8388608;
872	// redo with switch off
873	returnCode=gutsOfSolution ( givenDuals,givenPrimals,valuesPass);
874	if(primalTolerance_>testTolerance)
875	moreSpecialOptions_ \|= 8388608\|4194304;
876	if( !sumOfRelaxedPrimalInfeasibilities_)
877	doneFiddling=false; // over done it
878	}
879	}
880	}
881	}
882	}
883	return returnCode;
884	}
885	void
886	ClpSimplex::computePrimals ( const double * rowActivities,
887	const double * columnActivities)
888	{
889
890	//work space
891	CoinIndexedVector * workSpace = rowArray_[0];
892
893	CoinIndexedVector * arrayVector = rowArray_[1];
894	arrayVector->clear();
895	CoinIndexedVector * previousVector = rowArray_[2];
896	previousVector->clear();
897	// accumulate non basic stuff
898
899	int iRow;
900	// order is this way for scaling
901	if (columnActivities != columnActivityWork_)
902	ClpDisjointCopyN(columnActivities, numberColumns_, columnActivityWork_);
903	if (rowActivities != rowActivityWork_)
904	ClpDisjointCopyN(rowActivities, numberRows_, rowActivityWork_);
905	double * array = arrayVector->denseVector();
906	int * index = arrayVector->getIndices();
907	int number = 0;
908	const double * rhsOffset = matrix_->rhsOffset(this, false, true);
909	if (!rhsOffset) {
910	// Use whole matrix every time to make it easier for ClpMatrixBase
911	// So zero out basic
912	for (iRow = 0; iRow < numberRows_; iRow++) {
913	int iPivot = pivotVariable_[iRow];
914	assert (iPivot >= 0);
915	solution_[iPivot] = 0.0;
916	#ifdef CLP_INVESTIGATE
917	assert (getStatus(iPivot) == basic);
918	#endif
919	}
920	// Extended solution before "update"
921	matrix_->primalExpanded(this, 0);
922	times(-1.0, columnActivityWork_, array);
923	for (iRow = 0; iRow < numberRows_; iRow++) {
924	double value = array[iRow] + rowActivityWork_[iRow];
925	if (value) {
926	array[iRow] = value;
927	index[number++] = iRow;
928	} else {
929	array[iRow] = 0.0;
930	}
931	}
932	} else {
933	// we have an effective rhs lying around
934	// zero out basic (really just for slacks)
935	for (iRow = 0; iRow < numberRows_; iRow++) {
936	int iPivot = pivotVariable_[iRow];
937	solution_[iPivot] = 0.0;
938	}
939	for (iRow = 0; iRow < numberRows_; iRow++) {
940	double value = rhsOffset[iRow] + rowActivityWork_[iRow];
941	if (value) {
942	array[iRow] = value;
943	index[number++] = iRow;
944	} else {
945	array[iRow] = 0.0;
946	}
947	}
948	}
949	arrayVector->setNumElements(number);
950	#ifdef CLP_DEBUG
951	if (numberIterations_ == -3840) {
952	int i;
953	for (i = 0; i < numberRows_ + numberColumns_; i++)
954	printf("%d status %d\n", i, status_[i]);
955	printf("xxxxx1\n");
956	for (i = 0; i < numberRows_; i++)
957	if (array[i])
958	printf("%d rhs %g\n", i, array[i]);
959	printf("xxxxx2\n");
960	for (i = 0; i < numberRows_ + numberColumns_; i++)
961	if (getStatus(i) != basic)
962	printf("%d non basic %g %g %g\n", i, lower_[i], solution_[i], upper_[i]);
963	printf("xxxxx3\n");
964	}
965	#endif
966	// Ftran adjusted RHS and iterate to improve accuracy
967	double lastError = COIN_DBL_MAX;
968	int iRefine;
969	CoinIndexedVector * thisVector = arrayVector;
970	CoinIndexedVector * lastVector = previousVector;
971	#if 0
972	static double * xsave=NULL;
973	{
974	double * xx = thisVector->denseVector();
975	double largest=0.0;
976	int iLargest=-1;
977	for (int i=0;i<numberRows_;i++) {
978	if (fabs(xx[i])>largest) {
979	largest=fabs(xx[i]);
980	iLargest=i;
981	}
982	}
983	printf("largest incoming rhs %g on row %d\n",largest,iLargest);
984	}
985	if (numberIterations_<-40722) {
986	double * xx = thisVector->denseVector();
987	if (xsave) {
988	double * sol = xsave+numberRows_;
989	double largest=0.0;
990	int iLargest=-1;
991	for (int i=0;i<numberRows_;i++) {
992	if (fabs(xx[i]-xsave[i])>largest) {
993	largest=fabs(xx[i]-xsave[i]);
994	iLargest=i;
995	}
996	}
997	printf("error %g on row %d\n",largest,iLargest);
998	largest=0.0;
999	iLargest=-1;
1000	for (int i=0;i<numberColumns_;i++) {
1001	if (fabs(solution_[i]-sol[i])>largest) {
1002	largest=fabs(solution_[i]-sol[i]);
1003	iLargest=i;
1004	}
1005	}
1006	printf("error %g on col %d\n",largest,iLargest);
1007	} else {
1008	xsave=new double[2*numberRows_+numberColumns_];
1009	}
1010	memcpy(xsave,xx,numberRows_*sizeof(double));
1011	memcpy(xsave+numberRows_,solution_,(numberRows_+numberColumns_)*sizeof(double));
1012	}
1013	#endif
1014	//printf ("ZZ0 n before %d",number);
1015	if (number)
1016	factorization_->updateColumn(workSpace, thisVector);
1017	//printf(" - after %d\n",thisVector->getNumElements());
1018	double * work = workSpace->denseVector();
1019	#ifdef CLP_DEBUG
1020	if (numberIterations_ == -3840) {
1021	int i;
1022	for (i = 0; i < numberRows_; i++)
1023	if (array[i])
1024	printf("%d after rhs %g\n", i, array[i]);
1025	printf("xxxxx4\n");
1026	}
1027	#endif
1028	bool goodSolution = true;
1029	for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
1030
1031	int numberIn = thisVector->getNumElements();
1032	int * indexIn = thisVector->getIndices();
1033	double * arrayIn = thisVector->denseVector();
1034	// put solution in correct place
1035	if (!rhsOffset) {
1036	int j;
1037	for (j = 0; j < numberIn; j++) {
1038	iRow = indexIn[j];
1039	int iPivot = pivotVariable_[iRow];
1040	solution_[iPivot] = arrayIn[iRow];
1041	//assert (fabs(solution_[iPivot])<1.0e100);
1042	}
1043	} else {
1044	for (iRow = 0; iRow < numberRows_; iRow++) {
1045	int iPivot = pivotVariable_[iRow];
1046	solution_[iPivot] = arrayIn[iRow];
1047	//assert (fabs(solution_[iPivot])<1.0e100);
1048	}
1049	}
1050	// Extended solution after "update"
1051	matrix_->primalExpanded(this, 1);
1052	// check Ax == b (for all)
1053	// signal column generated matrix to just do basic (and gub)
1054	unsigned int saveOptions = specialOptions();
1055	setSpecialOptions(16);
1056	times(-1.0, columnActivityWork_, work);
1057	setSpecialOptions(saveOptions);
1058	largestPrimalError_ = 0.0;
1059	double multiplier = 131072.0;
1060	for (iRow = 0; iRow < numberRows_; iRow++) {
1061	double value = work[iRow] + rowActivityWork_[iRow];
1062	work[iRow] = value * multiplier;
1063	if (fabs(value) > largestPrimalError_) {
1064	largestPrimalError_ = fabs(value);
1065	}
1066	}
1067	if (largestPrimalError_ >= lastError) {
1068	// restore
1069	CoinIndexedVector * temp = thisVector;
1070	thisVector = lastVector;
1071	lastVector = temp;
1072	goodSolution = false;
1073	break;
1074	}
1075	if (iRefine < numberRefinements_ && largestPrimalError_ > 1.0e-10) {
1076	// try and make better
1077	// save this
1078	CoinIndexedVector * temp = thisVector;
1079	thisVector = lastVector;
1080	lastVector = temp;
1081	int * indexOut = thisVector->getIndices();
1082	int number = 0;
1083	array = thisVector->denseVector();
1084	thisVector->clear();
1085	for (iRow = 0; iRow < numberRows_; iRow++) {
1086	double value = work[iRow];
1087	if (value) {
1088	array[iRow] = value;
1089	indexOut[number++] = iRow;
1090	work[iRow] = 0.0;
1091	}
1092	}
1093	thisVector->setNumElements(number);
1094	lastError = largestPrimalError_;
1095	//printf ("ZZ%d n before %d",iRefine+1,number);
1096	factorization_->updateColumn(workSpace, thisVector);
1097	//printf(" - after %d\n",thisVector->getNumElements());
1098	multiplier = 1.0 / multiplier;
1099	double * previous = lastVector->denseVector();
1100	number = 0;
1101	for (iRow = 0; iRow < numberRows_; iRow++) {
1102	double value = previous[iRow] + multiplier * array[iRow];
1103	if (value) {
1104	array[iRow] = value;
1105	indexOut[number++] = iRow;
1106	} else {
1107	array[iRow] = 0.0;
1108	}
1109	}
1110	thisVector->setNumElements(number);
1111	} else {
1112	break;
1113	}
1114	}
1115
1116	// solution as accurate as we are going to get
1117	ClpFillN(work, numberRows_, 0.0);
1118	if (!goodSolution) {
1119	array = thisVector->denseVector();
1120	// put solution in correct place
1121	for (iRow = 0; iRow < numberRows_; iRow++) {
1122	int iPivot = pivotVariable_[iRow];
1123	solution_[iPivot] = array[iRow];
1124	//assert (fabs(solution_[iPivot])<1.0e100);
1125	}
1126	}
1127	arrayVector->clear();
1128	previousVector->clear();
1129	#ifdef CLP_DEBUG
1130	if (numberIterations_ == -3840) {
1131	exit(77);
1132	}
1133	#endif
1134	}
1135	// now dual side
1136	void
1137	ClpSimplex::computeDuals(double * givenDjs)
1138	{
1139	#ifndef SLIM_CLP
1140	if (objective_->type() == 1 \|\| !objective_->activated()) {
1141	#endif
1142	// Linear
1143	//work space
1144	CoinIndexedVector * workSpace = rowArray_[0];
1145
1146	CoinIndexedVector * arrayVector = rowArray_[1];
1147	arrayVector->clear();
1148	CoinIndexedVector * previousVector = rowArray_[2];
1149	previousVector->clear();
1150	int iRow;
1151	#ifdef CLP_DEBUG
1152	workSpace->checkClear();
1153	#endif
1154	double * array = arrayVector->denseVector();
1155	int * index = arrayVector->getIndices();
1156	int number = 0;
1157	if (!givenDjs) {
1158	for (iRow = 0; iRow < numberRows_; iRow++) {
1159	int iPivot = pivotVariable_[iRow];
1160	double value = cost_[iPivot];
1161	if (value) {
1162	array[iRow] = value;
1163	index[number++] = iRow;
1164	}
1165	}
1166	} else {
1167	// dual values pass - djs may not be zero
1168	for (iRow = 0; iRow < numberRows_; iRow++) {
1169	int iPivot = pivotVariable_[iRow];
1170	// make sure zero if done
1171	if (!pivoted(iPivot))
1172	givenDjs[iPivot] = 0.0;
1173	double value = cost_[iPivot] - givenDjs[iPivot];
1174	if (value) {
1175	array[iRow] = value;
1176	index[number++] = iRow;
1177	}
1178	}
1179	}
1180	arrayVector->setNumElements(number);
1181	// Extended duals before "updateTranspose"
1182	matrix_->dualExpanded(this, arrayVector, givenDjs, 0);
1183
1184	// Btran basic costs and get as accurate as possible
1185	double lastError = COIN_DBL_MAX;
1186	int iRefine;
1187	double * work = workSpace->denseVector();
1188	CoinIndexedVector * thisVector = arrayVector;
1189	CoinIndexedVector * lastVector = previousVector;
1190	factorization_->updateColumnTranspose(workSpace, thisVector);
1191
1192	for (iRefine = 0; iRefine < numberRefinements_ + 1; iRefine++) {
1193	// check basic reduced costs zero
1194	largestDualError_ = 0.0;
1195	if (!numberExtraRows_) {
1196	// Just basic
1197	int * index2 = workSpace->getIndices();
1198	// use reduced costs for slacks as work array
1199	double * work2 = reducedCostWork_ + numberColumns_;
1200	int numberStructurals = 0;
1201	for (iRow = 0; iRow < numberRows_; iRow++) {
1202	int iPivot = pivotVariable_[iRow];
1203	if (iPivot < numberColumns_)
1204	index2[numberStructurals++] = iPivot;
1205	}
1206	matrix_->listTransposeTimes(this, array, index2, numberStructurals, work2);
1207	numberStructurals = 0;
1208	if (!givenDjs) {
1209	for (iRow = 0; iRow < numberRows_; iRow++) {
1210	int iPivot = pivotVariable_[iRow];
1211	double value;
1212	if (iPivot >= numberColumns_) {
1213	// slack
1214	value = rowObjectiveWork_[iPivot-numberColumns_]
1215	+ array[iPivot-numberColumns_];
1216	} else {
1217	// column
1218	value = objectiveWork_[iPivot] - work2[numberStructurals++];
1219	}
1220	work[iRow] = value;
1221	if (fabs(value) > largestDualError_) {
1222	largestDualError_ = fabs(value);
1223	}
1224	}
1225	} else {
1226	for (iRow = 0; iRow < numberRows_; iRow++) {
1227	int iPivot = pivotVariable_[iRow];
1228	if (iPivot >= numberColumns_) {
1229	// slack
1230	work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
1231	+ array[iPivot-numberColumns_] - givenDjs[iPivot];
1232	} else {
1233	// column
1234	work[iRow] = objectiveWork_[iPivot] - work2[numberStructurals++]
1235	- givenDjs[iPivot];
1236	}
1237	if (fabs(work[iRow]) > largestDualError_) {
1238	largestDualError_ = fabs(work[iRow]);
1239	//assert (largestDualError_<1.0e-7);
1240	//if (largestDualError_>1.0e-7)
1241	//printf("large dual error %g\n",largestDualError_);
1242	}
1243	}
1244	}
1245	} else {
1246	// extra rows - be more careful
1247	#if 1
1248	// would be faster to do just for basic but this reduces code
1249	ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
1250	transposeTimes(-1.0, array, reducedCostWork_);
1251	#else
1252	// Just basic
1253	int * index2 = workSpace->getIndices();
1254	int numberStructurals = 0;
1255	for (iRow = 0; iRow < numberRows_; iRow++) {
1256	int iPivot = pivotVariable_[iRow];
1257	if (iPivot < numberColumns_)
1258	index2[numberStructurals++] = iPivot;
1259	}
1260	matrix_->listTransposeTimes(this, array, index2, numberStructurals, work);
1261	for (iRow = 0; iRow < numberStructurals; iRow++) {
1262	int iPivot = index2[iRow];
1263	reducedCostWork_[iPivot] = objectiveWork_[iPivot] - work[iRow];
1264	}
1265	#endif
1266	// update by duals on sets
1267	matrix_->dualExpanded(this, NULL, NULL, 1);
1268	if (!givenDjs) {
1269	for (iRow = 0; iRow < numberRows_; iRow++) {
1270	int iPivot = pivotVariable_[iRow];
1271	double value;
1272	if (iPivot >= numberColumns_) {
1273	// slack
1274	value = rowObjectiveWork_[iPivot-numberColumns_]
1275	+ array[iPivot-numberColumns_];
1276	} else {
1277	// column
1278	value = reducedCostWork_[iPivot];
1279	}
1280	work[iRow] = value;
1281	if (fabs(value) > largestDualError_) {
1282	largestDualError_ = fabs(value);
1283	}
1284	}
1285	} else {
1286	for (iRow = 0; iRow < numberRows_; iRow++) {
1287	int iPivot = pivotVariable_[iRow];
1288	if (iPivot >= numberColumns_) {
1289	// slack
1290	work[iRow] = rowObjectiveWork_[iPivot-numberColumns_]
1291	+ array[iPivot-numberColumns_] - givenDjs[iPivot];
1292	} else {
1293	// column
1294	work[iRow] = reducedCostWork_[iPivot] - givenDjs[iPivot];
1295	}
1296	if (fabs(work[iRow]) > largestDualError_) {
1297	largestDualError_ = fabs(work[iRow]);
1298	//assert (largestDualError_<1.0e-7);
1299	//if (largestDualError_>1.0e-7)
1300	//printf("large dual error %g\n",largestDualError_);
1301	}
1302	}
1303	}
1304	}
1305	if (largestDualError_ >= lastError) {
1306	// restore
1307	CoinIndexedVector * temp = thisVector;
1308	thisVector = lastVector;
1309	lastVector = temp;
1310	break;
1311	}
1312	if (iRefine < numberRefinements_ && largestDualError_ > 1.0e-10
1313	&& !givenDjs) {
1314	// try and make better
1315	// save this
1316	CoinIndexedVector * temp = thisVector;
1317	thisVector = lastVector;
1318	lastVector = temp;
1319	int * indexOut = thisVector->getIndices();
1320	int number = 0;
1321	array = thisVector->denseVector();
1322	thisVector->clear();
1323	double multiplier = 131072.0;
1324	for (iRow = 0; iRow < numberRows_; iRow++) {
1325	double value = multiplier * work[iRow];
1326	if (value) {
1327	array[iRow] = value;
1328	indexOut[number++] = iRow;
1329	work[iRow] = 0.0;
1330	}
1331	work[iRow] = 0.0;
1332	}
1333	thisVector->setNumElements(number);
1334	lastError = largestDualError_;
1335	factorization_->updateColumnTranspose(workSpace, thisVector);
1336	multiplier = 1.0 / multiplier;
1337	double * previous = lastVector->denseVector();
1338	number = 0;
1339	for (iRow = 0; iRow < numberRows_; iRow++) {
1340	double value = previous[iRow] + multiplier * array[iRow];
1341	if (value) {
1342	array[iRow] = value;
1343	indexOut[number++] = iRow;
1344	} else {
1345	array[iRow] = 0.0;
1346	}
1347	}
1348	thisVector->setNumElements(number);
1349	} else {
1350	break;
1351	}
1352	}
1353	// now look at dual solution
1354	array = thisVector->denseVector();
1355	for (iRow = 0; iRow < numberRows_; iRow++) {
1356	// slack
1357	double value = array[iRow];
1358	dual_[iRow] = value;
1359	value += rowObjectiveWork_[iRow];
1360	rowReducedCost_[iRow] = value;
1361	}
1362	// can use work if problem scaled (for better cache)
1363	ClpPackedMatrix* clpMatrix =
1364	dynamic_cast< ClpPackedMatrix*>(matrix_);
1365	double * saveRowScale = rowScale_;
1366	//double * saveColumnScale = columnScale_;
1367	if (scaledMatrix_) {
1368	rowScale_ = NULL;
1369	clpMatrix = scaledMatrix_;
1370	}
1371	if (clpMatrix && (clpMatrix->flags() & 2) == 0) {
1372	CoinIndexedVector * cVector = columnArray_[0];
1373	int * whichColumn = cVector->getIndices();
1374	assert (!cVector->getNumElements());
1375	int n = 0;
1376	for (int i = 0; i < numberColumns_; i++) {
1377	if (getColumnStatus(i) != basic) {
1378	whichColumn[n++] = i;
1379	reducedCostWork_[i] = objectiveWork_[i];
1380	} else {
1381	reducedCostWork_[i] = 0.0;
1382	}
1383	}
1384	if (numberRows_ > 4000)
1385	clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
1386	rowScale_, columnScale_, work);
1387	else
1388	clpMatrix->transposeTimesSubset(n, whichColumn, dual_, reducedCostWork_,
1389	rowScale_, columnScale_, NULL);
1390	} else {
1391	ClpDisjointCopyN(objectiveWork_, numberColumns_, reducedCostWork_);
1392	if (numberRows_ > 4000)
1393	matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
1394	rowScale_, columnScale_, work);
1395	else
1396	matrix_->transposeTimes(-1.0, dual_, reducedCostWork_,
1397	rowScale_, columnScale_, NULL);
1398	}
1399	rowScale_ = saveRowScale;
1400	//columnScale_ = saveColumnScale;
1401	ClpFillN(work, numberRows_, 0.0);
1402	// Extended duals and check dual infeasibility
1403	if (!matrix_->skipDualCheck() \|\| algorithm_ < 0 \|\| problemStatus_ != -2)
1404	matrix_->dualExpanded(this, NULL, NULL, 2);
1405	// If necessary - override results
1406	if (givenDjs) {
1407	// restore accurate duals
1408	CoinMemcpyN(dj_, (numberRows_ + numberColumns_), givenDjs);
1409	}
1410	arrayVector->clear();
1411	previousVector->clear();
1412	#ifndef SLIM_CLP
1413	} else {
1414	// Nonlinear
1415	objective_->reducedGradient(this, dj_, false);
1416	// get dual_ by moving from reduced costs for slacks
1417	CoinMemcpyN(dj_ + numberColumns_, numberRows_, dual_);
1418	}
1419	#endif
1420	}
1421	/* Given an existing factorization computes and checks
1422	primal and dual solutions. Uses input arrays for variables at
1423	bounds. Returns feasibility states */
1424	int ClpSimplex::getSolution ( const double * /rowActivities/,
1425	const double * /columnActivities/)
1426	{
1427	if (!factorization_->status()) {
1428	// put in standard form
1429	createRim(7 + 8 + 16 + 32, false, -1);
1430	if (pivotVariable_[0] < 0)
1431	internalFactorize(0);
1432	// do work
1433	gutsOfSolution ( NULL, NULL);
1434	// release extra memory
1435	deleteRim(0);
1436	}
1437	return factorization_->status();
1438	}
1439	/* Given an existing factorization computes and checks
1440	primal and dual solutions. Uses current problem arrays for
1441	bounds. Returns feasibility states */
1442	int ClpSimplex::getSolution ( )
1443	{
1444	double * rowActivities = new double[numberRows_];
1445	double * columnActivities = new double[numberColumns_];
1446	ClpDisjointCopyN ( rowActivityWork_, numberRows_ , rowActivities);
1447	ClpDisjointCopyN ( columnActivityWork_, numberColumns_ , columnActivities);
1448	int status = getSolution( rowActivities, columnActivities);
1449	delete [] rowActivities;
1450	delete [] columnActivities;
1451	return status;
1452	}
1453	// Factorizes using current basis. This is for external use
1454	// Return codes are as from ClpFactorization
1455	int ClpSimplex::factorize ()
1456	{
1457	// put in standard form
1458	createRim(7 + 8 + 16 + 32, false);
1459	// do work
1460	int status = internalFactorize(-1);
1461	// release extra memory
1462	deleteRim(0);
1463
1464	return status;
1465	}
1466	// Clean up status
1467	void
1468	ClpSimplex::cleanStatus()
1469	{
1470	int iRow, iColumn;
1471	int numberBasic = 0;
1472	// make row activities correct
1473	memset(rowActivityWork_, 0, numberRows_ * sizeof(double));
1474	times(1.0, columnActivityWork_, rowActivityWork_);
1475	if (!status_)
1476	createStatus();
1477	for (iRow = 0; iRow < numberRows_; iRow++) {
1478	if (getRowStatus(iRow) == basic)
1479	numberBasic++;
1480	else {
1481	setRowStatus(iRow, superBasic);
1482	// but put to bound if close
1483	if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
1484	<= primalTolerance_) {
1485	rowActivityWork_[iRow] = rowLowerWork_[iRow];
1486	setRowStatus(iRow, atLowerBound);
1487	} else if (fabs(rowActivityWork_[iRow] - rowUpperWork_[iRow])
1488	<= primalTolerance_) {
1489	rowActivityWork_[iRow] = rowUpperWork_[iRow];
1490	setRowStatus(iRow, atUpperBound);
1491	}
1492	}
1493	}
1494	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1495	if (getColumnStatus(iColumn) == basic) {
1496	if (numberBasic == numberRows_) {
1497	// take out of basis
1498	setColumnStatus(iColumn, superBasic);
1499	// but put to bound if close
1500	if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1501	<= primalTolerance_) {
1502	columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1503	setColumnStatus(iColumn, atLowerBound);
1504	} else if (fabs(columnActivityWork_[iColumn]
1505	- columnUpperWork_[iColumn])
1506	<= primalTolerance_) {
1507	columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1508	setColumnStatus(iColumn, atUpperBound);
1509	}
1510	} else
1511	numberBasic++;
1512	} else {
1513	setColumnStatus(iColumn, superBasic);
1514	// but put to bound if close
1515	if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1516	<= primalTolerance_) {
1517	columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1518	setColumnStatus(iColumn, atLowerBound);
1519	} else if (fabs(columnActivityWork_[iColumn]
1520	- columnUpperWork_[iColumn])
1521	<= primalTolerance_) {
1522	columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1523	setColumnStatus(iColumn, atUpperBound);
1524	}
1525	}
1526	}
1527	}
1528
1529	/* Factorizes using current basis.
1530	solveType - 1 iterating, 0 initial, -1 external
1531	- 2 then iterating but can throw out of basis
1532	If 10 added then in primal values pass
1533	Return codes are as from ClpFactorization unless initial factorization
1534	when total number of singularities is returned.
1535	Special case is numberRows_+1 -> all slack basis.
1536	*/
1537	int ClpSimplex::internalFactorize ( int solveType)
1538	{
1539	int iRow, iColumn;
1540	int totalSlacks = numberRows_;
1541	if (!status_)
1542	createStatus();
1543
1544	bool valuesPass = false;
1545	if (solveType >= 10) {
1546	valuesPass = true;
1547	solveType -= 10;
1548	}
1549	#ifdef CLP_DEBUG
1550	if (solveType > 0) {
1551	int numberFreeIn = 0, numberFreeOut = 0;
1552	double biggestDj = 0.0;
1553	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1554	switch(getColumnStatus(iColumn)) {
1555
1556	case basic:
1557	if (columnLower_[iColumn] < -largeValue_
1558	&& columnUpper_[iColumn] > largeValue_)
1559	numberFreeIn++;
1560	break;
1561	default:
1562	if (columnLower_[iColumn] < -largeValue_
1563	&& columnUpper_[iColumn] > largeValue_) {
1564	numberFreeOut++;
1565	biggestDj = CoinMax(fabs(dj_[iColumn]), biggestDj);
1566	}
1567	break;
1568	}
1569	}
1570	if (numberFreeIn + numberFreeOut)
1571	printf("%d in basis, %d out - largest dj %g\n",
1572	numberFreeIn, numberFreeOut, biggestDj);
1573	}
1574	#endif
1575	if (solveType <= 0) {
1576	// Make sure everything is clean
1577	for (iRow = 0; iRow < numberRows_; iRow++) {
1578	if(getRowStatus(iRow) == isFixed) {
1579	// double check fixed
1580	if (rowUpperWork_[iRow] > rowLowerWork_[iRow])
1581	setRowStatus(iRow, atLowerBound);
1582	} else if (getRowStatus(iRow) == isFree) {
1583	// may not be free after all
1584	if (rowLowerWork_[iRow] > -largeValue_ \|\| rowUpperWork_[iRow] < largeValue_)
1585	setRowStatus(iRow, superBasic);
1586	}
1587	}
1588	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1589	if(getColumnStatus(iColumn) == isFixed) {
1590	// double check fixed
1591	if (columnUpperWork_[iColumn] > columnLowerWork_[iColumn])
1592	setColumnStatus(iColumn, atLowerBound);
1593	} else if (getColumnStatus(iColumn) == isFree) {
1594	// may not be free after all
1595	if (columnLowerWork_[iColumn] > -largeValue_ \|\| columnUpperWork_[iColumn] < largeValue_)
1596	setColumnStatus(iColumn, superBasic);
1597	}
1598	}
1599	if (!valuesPass) {
1600	// not values pass so set to bounds
1601	bool allSlack = true;
1602	if (status_) {
1603	for (iRow = 0; iRow < numberRows_; iRow++) {
1604	if (getRowStatus(iRow) != basic) {
1605	allSlack = false;
1606	break;
1607	}
1608	}
1609	}
1610	if (!allSlack) {
1611	//#define CLP_INVESTIGATE2
1612	#ifdef CLP_INVESTIGATE3
1613	int numberTotal = numberRows_ + numberColumns_;
1614	double * saveSol = valuesPass ?
1615	CoinCopyOfArray(solution_, numberTotal) : NULL;
1616	#endif
1617	// set values from warm start (if sensible)
1618	int numberBasic = 0;
1619	for (iRow = 0; iRow < numberRows_; iRow++) {
1620	switch(getRowStatus(iRow)) {
1621
1622	case basic:
1623	numberBasic++;
1624	break;
1625	case atUpperBound:
1626	rowActivityWork_[iRow] = rowUpperWork_[iRow];
1627	if (rowActivityWork_[iRow] > largeValue_) {
1628	if (rowLowerWork_[iRow] > -largeValue_) {
1629	rowActivityWork_[iRow] = rowLowerWork_[iRow];
1630	setRowStatus(iRow, atLowerBound);
1631	} else {
1632	// say free
1633	setRowStatus(iRow, isFree);
1634	rowActivityWork_[iRow] = 0.0;
1635	}
1636	}
1637	break;
1638	case ClpSimplex::isFixed:
1639	case atLowerBound:
1640	rowActivityWork_[iRow] = rowLowerWork_[iRow];
1641	if (rowActivityWork_[iRow] < -largeValue_) {
1642	if (rowUpperWork_[iRow] < largeValue_) {
1643	rowActivityWork_[iRow] = rowUpperWork_[iRow];
1644	setRowStatus(iRow, atUpperBound);
1645	} else {
1646	// say free
1647	setRowStatus(iRow, isFree);
1648	rowActivityWork_[iRow] = 0.0;
1649	}
1650	}
1651	break;
1652	case isFree:
1653	break;
1654	// not really free - fall through to superbasic
1655	case superBasic:
1656	if (rowUpperWork_[iRow] > largeValue_) {
1657	if (rowLowerWork_[iRow] > -largeValue_) {
1658	rowActivityWork_[iRow] = rowLowerWork_[iRow];
1659	setRowStatus(iRow, atLowerBound);
1660	} else {
1661	// say free
1662	setRowStatus(iRow, isFree);
1663	rowActivityWork_[iRow] = 0.0;
1664	}
1665	} else {
1666	if (rowLowerWork_[iRow] > -largeValue_) {
1667	// set to nearest
1668	if (fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])
1669	< fabs(rowActivityWork_[iRow] - rowLowerWork_[iRow])) {
1670	rowActivityWork_[iRow] = rowLowerWork_[iRow];
1671	setRowStatus(iRow, atLowerBound);
1672	} else {
1673	rowActivityWork_[iRow] = rowUpperWork_[iRow];
1674	setRowStatus(iRow, atUpperBound);
1675	}
1676	} else {
1677	rowActivityWork_[iRow] = rowUpperWork_[iRow];
1678	setRowStatus(iRow, atUpperBound);
1679	}
1680	}
1681	break;
1682	}
1683	}
1684	totalSlacks = numberBasic;
1685
1686	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1687	switch(getColumnStatus(iColumn)) {
1688
1689	case basic:
1690	if (numberBasic == maximumBasic_) {
1691	// take out of basis
1692	if (columnLowerWork_[iColumn] > -largeValue_) {
1693	if (columnActivityWork_[iColumn] - columnLowerWork_[iColumn] <
1694	columnUpperWork_[iColumn] - columnActivityWork_[iColumn]) {
1695	columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1696	setColumnStatus(iColumn, atLowerBound);
1697	} else {
1698	columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1699	setColumnStatus(iColumn, atUpperBound);
1700	}
1701	} else if (columnUpperWork_[iColumn] < largeValue_) {
1702	columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1703	setColumnStatus(iColumn, atUpperBound);
1704	} else {
1705	columnActivityWork_[iColumn] = 0.0;
1706	setColumnStatus(iColumn, isFree);
1707	}
1708	} else {
1709	numberBasic++;
1710	}
1711	break;
1712	case atUpperBound:
1713	columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1714	if (columnActivityWork_[iColumn] > largeValue_) {
1715	if (columnLowerWork_[iColumn] < -largeValue_) {
1716	columnActivityWork_[iColumn] = 0.0;
1717	setColumnStatus(iColumn, isFree);
1718	} else {
1719	columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1720	setColumnStatus(iColumn, atLowerBound);
1721	}
1722	}
1723	break;
1724	case isFixed:
1725	case atLowerBound:
1726	columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1727	if (columnActivityWork_[iColumn] < -largeValue_) {
1728	if (columnUpperWork_[iColumn] > largeValue_) {
1729	columnActivityWork_[iColumn] = 0.0;
1730	setColumnStatus(iColumn, isFree);
1731	} else {
1732	columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1733	setColumnStatus(iColumn, atUpperBound);
1734	}
1735	}
1736	break;
1737	case isFree:
1738	break;
1739	// not really free - fall through to superbasic
1740	case superBasic:
1741	if (columnUpperWork_[iColumn] > largeValue_) {
1742	if (columnLowerWork_[iColumn] > -largeValue_) {
1743	columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1744	setColumnStatus(iColumn, atLowerBound);
1745	} else {
1746	// say free
1747	setColumnStatus(iColumn, isFree);
1748	columnActivityWork_[iColumn] = 0.0;
1749	}
1750	} else {
1751	if (columnLowerWork_[iColumn] > -largeValue_) {
1752	// set to nearest
1753	if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1754	< fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])) {
1755	columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1756	setColumnStatus(iColumn, atLowerBound);
1757	} else {
1758	columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1759	setColumnStatus(iColumn, atUpperBound);
1760	}
1761	} else {
1762	columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1763	setColumnStatus(iColumn, atUpperBound);
1764	}
1765	}
1766	break;
1767	}
1768	}
1769	#ifdef CLP_INVESTIGATE3
1770	if (saveSol) {
1771	int numberChanged = 0;
1772	double largestChanged = 0.0;
1773	for (int i = 0; i < numberTotal; i++) {
1774	double difference = fabs(solution_[i] - saveSol[i]);
1775	if (difference > 1.0e-7) {
1776	numberChanged++;
1777	if (difference > largestChanged)
1778	largestChanged = difference;
1779	}
1780	}
1781	if (numberChanged)
1782	printf("%d changed, largest %g\n", numberChanged, largestChanged);
1783	delete [] saveSol;
1784	}
1785	#endif
1786	#if 0
1787	if (numberBasic < numberRows_) {
1788	// add some slacks in case odd warmstart
1789	#ifdef CLP_INVESTIGATE
1790	printf("BAD %d basic, %d rows %d slacks\n",
1791	numberBasic, numberRows_, totalSlacks);
1792	#endif
1793	int iRow = numberRows_ - 1;
1794	while (numberBasic < numberRows_) {
1795	if (getRowStatus(iRow) != basic) {
1796	setRowStatus(iRow, basic);
1797	numberBasic++;
1798	totalSlacks++;
1799	iRow--;
1800	} else {
1801	break;
1802	}
1803	}
1804	}
1805	#endif
1806	} else {
1807	// all slack basis
1808	int numberBasic = 0;
1809	if (!status_) {
1810	createStatus();
1811	}
1812	for (iRow = 0; iRow < numberRows_; iRow++) {
1813	double lower = rowLowerWork_[iRow];
1814	double upper = rowUpperWork_[iRow];
1815	if (lower > -largeValue_ \|\| upper < largeValue_) {
1816	if (fabs(lower) <= fabs(upper)) {
1817	rowActivityWork_[iRow] = lower;
1818	} else {
1819	rowActivityWork_[iRow] = upper;
1820	}
1821	} else {
1822	rowActivityWork_[iRow] = 0.0;
1823	}
1824	setRowStatus(iRow, basic);
1825	numberBasic++;
1826	}
1827	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1828	double lower = columnLowerWork_[iColumn];
1829	double upper = columnUpperWork_[iColumn];
1830	double big_bound = largeValue_;
1831	if (lower > -big_bound \|\| upper < big_bound) {
1832	if ((getColumnStatus(iColumn) == atLowerBound &&
1833	columnActivityWork_[iColumn] == lower) \|\|
1834	(getColumnStatus(iColumn) == atUpperBound &&
1835	columnActivityWork_[iColumn] == upper)) {
1836	// status looks plausible
1837	} else {
1838	// set to sensible
1839	if (fabs(lower) <= fabs(upper)) {
1840	setColumnStatus(iColumn, atLowerBound);
1841	columnActivityWork_[iColumn] = lower;
1842	} else {
1843	setColumnStatus(iColumn, atUpperBound);
1844	columnActivityWork_[iColumn] = upper;
1845	}
1846	}
1847	} else {
1848	setColumnStatus(iColumn, isFree);
1849	columnActivityWork_[iColumn] = 0.0;
1850	}
1851	}
1852	}
1853	} else {
1854	// values pass has less coding
1855	// make row activities correct and clean basis a bit
1856	cleanStatus();
1857	if (status_) {
1858	int numberBasic = 0;
1859	for (iRow = 0; iRow < numberRows_; iRow++) {
1860	if (getRowStatus(iRow) == basic)
1861	numberBasic++;
1862	}
1863	totalSlacks = numberBasic;
1864	#if 0
1865	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1866	if (getColumnStatus(iColumn) == basic)
1867	numberBasic++;
1868	}
1869	#endif
1870	} else {
1871	// all slack basis
1872	int numberBasic = 0;
1873	if (!status_) {
1874	createStatus();
1875	}
1876	for (iRow = 0; iRow < numberRows_; iRow++) {
1877	setRowStatus(iRow, basic);
1878	numberBasic++;
1879	}
1880	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1881	setColumnStatus(iColumn, superBasic);
1882	// but put to bound if close
1883	if (fabs(columnActivityWork_[iColumn] - columnLowerWork_[iColumn])
1884	<= primalTolerance_) {
1885	columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1886	setColumnStatus(iColumn, atLowerBound);
1887	} else if (fabs(columnActivityWork_[iColumn]
1888	- columnUpperWork_[iColumn])
1889	<= primalTolerance_) {
1890	columnActivityWork_[iColumn] = columnUpperWork_[iColumn];
1891	setColumnStatus(iColumn, atUpperBound);
1892	}
1893	}
1894	}
1895	}
1896	numberRefinements_ = 1;
1897	// set fixed if they are
1898	for (iRow = 0; iRow < numberRows_; iRow++) {
1899	if (getRowStatus(iRow) != basic ) {
1900	if (rowLowerWork_[iRow] == rowUpperWork_[iRow]) {
1901	rowActivityWork_[iRow] = rowLowerWork_[iRow];
1902	setRowStatus(iRow, isFixed);
1903	}
1904	}
1905	}
1906	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1907	if (getColumnStatus(iColumn) != basic ) {
1908	if (columnLowerWork_[iColumn] == columnUpperWork_[iColumn]) {
1909	columnActivityWork_[iColumn] = columnLowerWork_[iColumn];
1910	setColumnStatus(iColumn, isFixed);
1911	}
1912	}
1913	}
1914	}
1915	//for (iRow=0;iRow<numberRows_+numberColumns_;iRow++) {
1916	//if (fabs(solution_[iRow])>1.0e10) {
1917	// printf("large %g at %d - status %d\n",
1918	// solution_[iRow],iRow,status_[iRow]);
1919	//}
1920	//}
1921	# if 0 //ndef _MSC_VER
1922	// The local static var k is a problem when trying to build a DLL. Since this is
1923	// just for debugging (likely done on *nix), just hide it from Windows
1924	// -- lh, 101016 --
1925	if (0) {
1926	static int k = 0;
1927	printf("start basis\n");
1928	int i;
1929	for (i = 0; i < numberRows_; i++)
1930	printf ("xx %d %d\n", i, pivotVariable_[i]);
1931	for (i = 0; i < numberRows_ + numberColumns_; i++)
1932	if (getColumnStatus(i) == basic)
1933	printf ("yy %d basic\n", i);
1934	if (k > 20)
1935	exit(0);
1936	k++;
1937	}
1938	# endif
1939	#if 0 //ndef NDEBUG
1940	// Make sure everything is clean
1941	double sumOutside=0.0;
1942	int numberOutside=0;
1943	//double sumOutsideLarge=0.0;
1944	int numberOutsideLarge=0;
1945	double sumInside=0.0;
1946	int numberInside=0;
1947	//double sumInsideLarge=0.0;
1948	int numberInsideLarge=0;
1949	int numberTotal=numberRows_+numberColumns_;
1950	for (int iSequence = 0; iSequence < numberTotal; iSequence++) {
1951	if(getStatus(iSequence) == isFixed) {
1952	// double check fixed
1953	assert (upper_[iSequence] == lower_[iSequence]);
1954	assert (fabs(solution_[iSequence]-lower_[iSequence])<primalTolerance_);
1955	} else if (getStatus(iSequence) == isFree) {
1956	assert (upper_[iSequence] == COIN_DBL_MAX && lower_[iSequence]==-COIN_DBL_MAX);
1957	} else if (getStatus(iSequence) == atLowerBound) {
1958	assert (fabs(solution_[iSequence]-lower_[iSequence])<1000.0*primalTolerance_);
1959	if (solution_[iSequence]<lower_[iSequence]) {
1960	numberOutside++;
1961	sumOutside-=solution_[iSequence]-lower_[iSequence];
1962	if (solution_[iSequence]<lower_[iSequence]-primalTolerance_)
1963	numberOutsideLarge++;
1964	} else if (solution_[iSequence]>lower_[iSequence]) {
1965	numberInside++;
1966	sumInside+=solution_[iSequence]-lower_[iSequence];
1967	if (solution_[iSequence]>lower_[iSequence]+primalTolerance_)
1968	numberInsideLarge++;
1969	}
1970	} else if (getStatus(iSequence) == atUpperBound) {
1971	assert (fabs(solution_[iSequence]-upper_[iSequence])<1000.0*primalTolerance_);
1972	if (solution_[iSequence]>upper_[iSequence]) {
1973	numberOutside++;
1974	sumOutside+=solution_[iSequence]-upper_[iSequence];
1975	if (solution_[iSequence]>upper_[iSequence]+primalTolerance_)
1976	numberOutsideLarge++;
1977	} else if (solution_[iSequence]<upper_[iSequence]) {
1978	numberInside++;
1979	sumInside-=solution_[iSequence]-upper_[iSequence];
1980	if (solution_[iSequence]<upper_[iSequence]-primalTolerance_)
1981	numberInsideLarge++;
1982	}
1983	} else if (getStatus(iSequence) == superBasic) {
1984	//assert (!valuesPass);
1985	}
1986	}
1987	if (numberInside+numberOutside)
1988	printf("%d outside summing to %g (%d large), %d inside summing to %g (%d large)\n",
1989	numberOutside,sumOutside,numberOutsideLarge,
1990	numberInside,sumInside,numberInsideLarge);
1991	#endif
1992	int status = factorization_->factorize(this, solveType, valuesPass);
1993	if (status) {
1994	handler_->message(CLP_SIMPLEX_BADFACTOR, messages_)
1995	<< status
1996	<< CoinMessageEol;
1997	#ifdef CLP_USEFUL_PRINTOUT
1998	printf("Basis singular - pivot tolerance %g\n",
1999	factorization_->pivotTolerance());
2000	#endif
2001	return -1;
2002	} else if (!solveType) {
2003	// Initial basis - return number of singularities
2004	int numberSlacks = 0;
2005	for (iRow = 0; iRow < numberRows_; iRow++) {
2006	if (getRowStatus(iRow) == basic)
2007	numberSlacks++;
2008	}
2009	status = CoinMax(numberSlacks - totalSlacks, 0);
2010	// special case if all slack
2011	if (numberSlacks == numberRows_) {
2012	status = numberRows_ + 1;
2013	}
2014	}
2015
2016	// sparse methods
2017	//if (factorization_->sparseThreshold()) {
2018	// get default value
2019	factorization_->sparseThreshold(0);
2020	if (!(moreSpecialOptions_&1024))
2021	factorization_->goSparse();
2022	//}
2023
2024	return status;
2025	}
2026	/*
2027	This does basis housekeeping and does values for in/out variables.
2028	Can also decide to re-factorize
2029	*/
2030	int
2031	ClpSimplex::housekeeping(double objectiveChange)
2032	{
2033	// save value of incoming and outgoing
2034	double oldIn = solution_[sequenceIn_];
2035	double oldOut = solution_[sequenceOut_];
2036	numberIterations_++;
2037	changeMade_++; // something has happened
2038	// incoming variable
2039	if (handler_->logLevel() > 7) {
2040	//if (handler_->detail(CLP_SIMPLEX_HOUSE1,messages_)<100) {
2041	handler_->message(CLP_SIMPLEX_HOUSE1, messages_)
2042	<< directionOut_
2043	<< directionIn_ << theta_
2044	<< dualOut_ << dualIn_ << alpha_
2045	<< CoinMessageEol;
2046	if (getStatus(sequenceIn_) == isFree) {
2047	handler_->message(CLP_SIMPLEX_FREEIN, messages_)
2048	<< sequenceIn_
2049	<< CoinMessageEol;
2050	}
2051	}
2052	#if 0
2053	printf("h1 %d %d %g %g %g %g",
2054	directionOut_
2055	, directionIn_, theta_
2056	, dualOut_, dualIn_, alpha_);
2057	#endif
2058	// change of incoming
2059	char rowcol[] = {'R', 'C'};
2060	if (pivotRow_ >= 0)
2061	pivotVariable_[pivotRow_] = sequenceIn();
2062	if (upper_[sequenceIn_] > 1.0e20 && lower_[sequenceIn_] < -1.0e20)
2063	progressFlag_ \|= 2; // making real progress
2064	solution_[sequenceIn_] = valueIn_;
2065	if (upper_[sequenceOut_] - lower_[sequenceOut_] < 1.0e-12)
2066	progressFlag_ \|= 1; // making real progress
2067	if (sequenceIn_ != sequenceOut_) {
2068	if (alphaAccuracy_ > 0.0) {
2069	double value = fabs(alpha_);
2070	if (value > 1.0)
2071	alphaAccuracy_ *= value;
2072	else
2073	alphaAccuracy_ /= value;
2074	}
2075	//assert( getStatus(sequenceOut_)== basic);
2076	setStatus(sequenceIn_, basic);
2077	if (upper_[sequenceOut_] - lower_[sequenceOut_] > 0) {
2078	// As Nonlinear costs may have moved bounds (to more feasible)
2079	// Redo using value
2080	if (fabs(valueOut_ - lower_[sequenceOut_]) < fabs(valueOut_ - upper_[sequenceOut_])) {
2081	// going to lower
2082	setStatus(sequenceOut_, atLowerBound);
2083	oldOut = lower_[sequenceOut_];
2084	} else {
2085	// going to upper
2086	setStatus(sequenceOut_, atUpperBound);
2087	oldOut = upper_[sequenceOut_];
2088	}
2089	} else {
2090	// fixed
2091	setStatus(sequenceOut_, isFixed);
2092	}
2093	solution_[sequenceOut_] = valueOut_;
2094	} else {
2095	//if (objective_->type()<2)
2096	//assert (fabs(theta_)>1.0e-13);
2097	// flip from bound to bound
2098	// As Nonlinear costs may have moved bounds (to more feasible)
2099	// Redo using value
2100	if (fabs(valueIn_ - lower_[sequenceIn_]) < fabs(valueIn_ - upper_[sequenceIn_])) {
2101	// as if from upper bound
2102	setStatus(sequenceIn_, atLowerBound);
2103	} else {
2104	// as if from lower bound
2105	setStatus(sequenceIn_, atUpperBound);
2106	}
2107	}
2108
2109	// Update hidden stuff e.g. effective RHS and gub
2110	int invertNow=matrix_->updatePivot(this, oldIn, oldOut);
2111	objectiveValue_ += objectiveChange / (objectiveScale_ * rhsScale_);
2112	if (handler_->logLevel() > 7) {
2113	//if (handler_->detail(CLP_SIMPLEX_HOUSE2,messages_)<100) {
2114	handler_->message(CLP_SIMPLEX_HOUSE2, messages_)
2115	<< numberIterations_ << objectiveValue()
2116	<< rowcol[isColumn(sequenceIn_)] << sequenceWithin(sequenceIn_)
2117	<< rowcol[isColumn(sequenceOut_)] << sequenceWithin(sequenceOut_);
2118	handler_->printing(algorithm_ < 0) << dualOut_ << theta_;
2119	handler_->printing(algorithm_ > 0) << dualIn_ << theta_;
2120	handler_->message() << CoinMessageEol;
2121	}
2122	#if 0
2123	if (numberIterations_ > 10000)
2124	printf(" it %d %g %c%d %c%d\n"
2125	, numberIterations_, objectiveValue()
2126	, rowcol[isColumn(sequenceIn_)], sequenceWithin(sequenceIn_)
2127	, rowcol[isColumn(sequenceOut_)], sequenceWithin(sequenceOut_));
2128	#endif
2129	if (trustedUserPointer_ && trustedUserPointer_->typeStruct == 1) {
2130	if (algorithm_ > 0 && integerType_ && !nonLinearCost_->numberInfeasibilities()) {
2131	if (fabs(theta_) > 1.0e-6 \|\| !numberIterations_) {
2132	// For saving solutions
2133	typedef struct {
2134	int numberSolutions;
2135	int maximumSolutions;
2136	int numberColumns;
2137	double ** solution;
2138	int * numberUnsatisfied;
2139	} clpSolution;
2140	clpSolution * solution = reinterpret_cast<clpSolution *> (trustedUserPointer_->data);
2141	if (solution->numberSolutions == solution->maximumSolutions) {
2142	int n = solution->maximumSolutions;
2143	int n2 = (n * 3) / 2 + 10;
2144	solution->maximumSolutions = n2;
2145	double ** temp = new double * [n2];
2146	for (int i = 0; i < n; i++)
2147	temp[i] = solution->solution[i];
2148	delete [] solution->solution;
2149	solution->solution = temp;
2150	int * tempN = new int [n2];
2151	for (int i = 0; i < n; i++)
2152	tempN[i] = solution->numberUnsatisfied[i];
2153	delete [] solution->numberUnsatisfied;
2154	solution->numberUnsatisfied = tempN;
2155	}
2156	assert (numberColumns_ == solution->numberColumns);
2157	double * sol = new double [numberColumns_];
2158	solution->solution[solution->numberSolutions] = sol;
2159	int numberFixed = 0;
2160	int numberUnsat = 0;
2161	int numberSat = 0;
2162	double sumUnsat = 0.0;
2163	double tolerance = 10.0 * primalTolerance_;
2164	double mostAway = 0.0;
2165	for (int i = 0; i < numberColumns_; i++) {
2166	// Save anyway
2167	sol[i] = columnScale_ ? solution_[i] * columnScale_[i] : solution_[i];
2168	// rest is optional
2169	if (upper_[i] > lower_[i]) {
2170	double value = solution_[i];
2171	if (value > lower_[i] + tolerance &&
2172	value < upper_[i] - tolerance && integerType_[i]) {
2173	// may have to modify value if scaled
2174	if (columnScale_)
2175	value *= columnScale_[i];
2176	double closest = floor(value + 0.5);
2177	// problem may be perturbed so relax test
2178	if (fabs(value - closest) > 1.0e-4) {
2179	numberUnsat++;
2180	sumUnsat += fabs(value - closest);
2181	if (mostAway < fabs(value - closest)) {
2182	mostAway = fabs(value - closest);
2183	}
2184	} else {
2185	numberSat++;
2186	}
2187	} else {
2188	numberSat++;
2189	}
2190	} else {
2191	numberFixed++;
2192	}
2193	}
2194	solution->numberUnsatisfied[solution->numberSolutions++] = numberUnsat;
2195	COIN_DETAIL_PRINT(printf("iteration %d, %d unsatisfied (%g,%g), %d fixed, %d satisfied\n",
2196	numberIterations_, numberUnsat, sumUnsat, mostAway, numberFixed, numberSat));
2197	}
2198	}
2199	}
2200	if (hitMaximumIterations())
2201	return 2;
2202	#if 1
2203	//if (numberIterations_>14000)
2204	//handler_->setLogLevel(63);
2205	//if (numberIterations_>24000)
2206	//exit(77);
2207	// check for small cycles
2208	int in = sequenceIn_;
2209	int out = sequenceOut_;
2210	matrix_->correctSequence(this, in, out);
2211	int cycle = progress_.cycle(in, out,
2212	directionIn_, directionOut_);
2213	if (cycle > 0 && objective_->type() < 2 && matrix_->type() < 15) {
2214	//if (cycle>0) {
2215	if (handler_->logLevel() >= 63)
2216	printf("Cycle of %d\n", cycle);
2217	// reset
2218	progress_.startCheck();
2219	double random = randomNumberGenerator_.randomDouble();
2220	int extra = static_cast<int> (9.999 * random);
2221	int off[] = {1, 1, 1, 1, 2, 2, 2, 3, 3, 4};
2222	if (factorization_->pivots() > cycle) {
2223	forceFactorization_ = CoinMax(1, cycle - off[extra]);
2224	} else {
2225	/* need to reject something
2226	should be better if don't reject incoming
2227	as it is in basis */
2228	int iSequence;
2229	//if (algorithm_ > 0)
2230	// iSequence = sequenceIn_;
2231	//else
2232	iSequence = sequenceOut_;
2233	char x = isColumn(iSequence) ? 'C' : 'R';
2234	if (handler_->logLevel() >= 63)
2235	handler_->message(CLP_SIMPLEX_FLAG, messages_)
2236	<< x << sequenceWithin(iSequence)
2237	<< CoinMessageEol;
2238	setFlagged(iSequence);
2239	//printf("flagging %d\n",iSequence);
2240	}
2241	return 1;
2242	}
2243	#endif
2244	// only time to re-factorize if one before real time
2245	// this is so user won't be surprised that maximumPivots has exact meaning
2246	int numberPivots = factorization_->pivots();
2247	int maximumPivots = factorization_->maximumPivots();
2248	int numberDense = factorization_->numberDense();
2249	bool dontInvert = ((specialOptions_ & 16384) != 0 && numberIterations_ * 3 >
2250	2 * maximumIterations());
2251	if (numberPivots == maximumPivots \|\|
2252	maximumPivots < 2) {
2253	// If dense then increase
2254	if (maximumPivots > 100 && numberDense > 1.5 * maximumPivots
2255	&& false) {
2256	factorization_->maximumPivots(numberDense);
2257	dualRowPivot_->maximumPivotsChanged();
2258	primalColumnPivot_->maximumPivotsChanged();
2259	// and redo arrays
2260	for (int iRow = 0; iRow < 4; iRow++) {
2261	int length = rowArray_[iRow]->capacity() + numberDense - maximumPivots;
2262	rowArray_[iRow]->reserve(length);
2263	}
2264	}
2265	#if CLP_FACTORIZATION_NEW_TIMING>1
2266	factorization_->statsRefactor('M');
2267	#endif
2268	return 1;
2269	} else if ((factorization_->timeToRefactorize() && !dontInvert)
2270	\|\|invertNow) {
2271	//printf("ret after %d pivots\n",factorization_->pivots());
2272	#if CLP_FACTORIZATION_NEW_TIMING>1
2273	factorization_->statsRefactor('T');
2274	#endif
2275	return 1;
2276	} else if (forceFactorization_ > 0 &&
2277	factorization_->pivots() == forceFactorization_) {
2278	// relax
2279	forceFactorization_ = (3 + 5 * forceFactorization_) / 4;
2280	if (forceFactorization_ > factorization_->maximumPivots())
2281	forceFactorization_ = -1; //off
2282	#if CLP_FACTORIZATION_NEW_TIMING>1
2283	factorization_->statsRefactor('F');
2284	#endif
2285	return 1;
2286	} else if (numberIterations_ > 1000 + 10 * (numberRows_ + (numberColumns_ >> 2)) && matrix_->type()<15) {
2287	double random = randomNumberGenerator_.randomDouble();
2288	while (random<0.45)
2289	random *= 2.0;
2290	int maxNumber = (forceFactorization_ < 0) ? maximumPivots : CoinMin(forceFactorization_, maximumPivots);
2291	if (factorization_->pivots() >= random * maxNumber) {
2292	return 1;
2293	} else if (numberIterations_ > 1000000 + 10 * (numberRows_ + (numberColumns_ >> 2)) &&
2294	numberIterations_ < 1001000 + 10 * (numberRows_ + (numberColumns_ >> 2))) {
2295	return 1;
2296	} else {
2297	// carry on iterating
2298	return 0;
2299	}
2300	} else {
2301	// carry on iterating
2302	return 0;
2303	}
2304	}
2305	// Copy constructor.
2306	ClpSimplex::ClpSimplex(const ClpSimplex &rhs, int scalingMode) :
2307	ClpModel(rhs, scalingMode),
2308	bestPossibleImprovement_(0.0),
2309	zeroTolerance_(1.0e-13),
2310	columnPrimalSequence_(-2),
2311	rowPrimalSequence_(-2),
2312	bestObjectiveValue_(rhs.bestObjectiveValue_),
2313	moreSpecialOptions_(2),
2314	baseIteration_(0),
2315	primalToleranceToGetOptimal_(-1.0),
2316	largeValue_(1.0e15),
2317	largestPrimalError_(0.0),
2318	largestDualError_(0.0),
2319	alphaAccuracy_(-1.0),
2320	dualBound_(1.0e10),
2321	alpha_(0.0),
2322	theta_(0.0),
2323	lowerIn_(0.0),
2324	valueIn_(0.0),
2325	upperIn_(-COIN_DBL_MAX),
2326	dualIn_(0.0),
2327	lowerOut_(-1),
2328	valueOut_(-1),
2329	upperOut_(-1),
2330	dualOut_(-1),
2331	dualTolerance_(1.0e-7),
2332	primalTolerance_(1.0e-7),
2333	sumDualInfeasibilities_(0.0),
2334	sumPrimalInfeasibilities_(0.0),
2335	infeasibilityCost_(1.0e10),
2336	sumOfRelaxedDualInfeasibilities_(0.0),
2337	sumOfRelaxedPrimalInfeasibilities_(0.0),
2338	acceptablePivot_(1.0e-8),
2339	lower_(NULL),
2340	rowLowerWork_(NULL),
2341	columnLowerWork_(NULL),
2342	upper_(NULL),
2343	rowUpperWork_(NULL),
2344	columnUpperWork_(NULL),
2345	cost_(NULL),
2346	rowObjectiveWork_(NULL),
2347	objectiveWork_(NULL),
2348	sequenceIn_(-1),
2349	directionIn_(-1),
2350	sequenceOut_(-1),
2351	directionOut_(-1),
2352	pivotRow_(-1),
2353	lastGoodIteration_(-100),
2354	dj_(NULL),
2355	rowReducedCost_(NULL),
2356	reducedCostWork_(NULL),
2357	solution_(NULL),
2358	rowActivityWork_(NULL),
2359	columnActivityWork_(NULL),
2360	numberDualInfeasibilities_(0),
2361	numberDualInfeasibilitiesWithoutFree_(0),
2362	numberPrimalInfeasibilities_(100),
2363	numberRefinements_(0),
2364	pivotVariable_(NULL),
2365	factorization_(NULL),
2366	savedSolution_(NULL),
2367	numberTimesOptimal_(0),
2368	disasterArea_(NULL),
2369	changeMade_(1),
2370	algorithm_(0),
2371	forceFactorization_(-1),
2372	perturbation_(100),
2373	nonLinearCost_(NULL),
2374	lastBadIteration_(-999999),
2375	lastFlaggedIteration_(-999999),
2376	numberFake_(0),
2377	numberChanged_(0),
2378	progressFlag_(0),
2379	firstFree_(-1),
2380	numberExtraRows_(0),
2381	maximumBasic_(0),
2382	dontFactorizePivots_(0),
2383	incomingInfeasibility_(1.0),
2384	allowedInfeasibility_(10.0),
2385	automaticScale_(0),
2386	maximumPerturbationSize_(0),
2387	perturbationArray_(NULL),
2388	baseModel_(NULL)
2389	#ifdef ABC_INHERIT
2390	,abcSimplex_(NULL),
2391	abcState_(0)
2392	#endif
2393	{
2394	int i;
2395	for (i = 0; i < 6; i++) {
2396	rowArray_[i] = NULL;
2397	columnArray_[i] = NULL;
2398	}
2399	for (i = 0; i < 4; i++) {
2400	spareIntArray_[i] = 0;
2401	spareDoubleArray_[i] = 0.0;
2402	}
2403	saveStatus_ = NULL;
2404	factorization_ = NULL;
2405	dualRowPivot_ = NULL;
2406	primalColumnPivot_ = NULL;
2407	gutsOfDelete(0);
2408	delete nonLinearCost_;
2409	nonLinearCost_ = NULL;
2410	gutsOfCopy(rhs);
2411	solveType_ = 1; // say simplex based life form
2412	}
2413	// Copy constructor from model
2414	ClpSimplex::ClpSimplex(const ClpModel &rhs, int scalingMode) :
2415	ClpModel(rhs, scalingMode),
2416	bestPossibleImprovement_(0.0),
2417	zeroTolerance_(1.0e-13),
2418	columnPrimalSequence_(-2),
2419	rowPrimalSequence_(-2),
2420	bestObjectiveValue_(-COIN_DBL_MAX),
2421	moreSpecialOptions_(2),
2422	baseIteration_(0),
2423	primalToleranceToGetOptimal_(-1.0),
2424	largeValue_(1.0e15),
2425	largestPrimalError_(0.0),
2426	largestDualError_(0.0),
2427	alphaAccuracy_(-1.0),
2428	dualBound_(1.0e10),
2429	alpha_(0.0),
2430	theta_(0.0),
2431	lowerIn_(0.0),
2432	valueIn_(0.0),
2433	upperIn_(-COIN_DBL_MAX),
2434	dualIn_(0.0),
2435	lowerOut_(-1),
2436	valueOut_(-1),
2437	upperOut_(-1),
2438	dualOut_(-1),
2439	dualTolerance_(1.0e-7),
2440	primalTolerance_(1.0e-7),
2441	sumDualInfeasibilities_(0.0),
2442	sumPrimalInfeasibilities_(0.0),
2443	infeasibilityCost_(1.0e10),
2444	sumOfRelaxedDualInfeasibilities_(0.0),
2445	sumOfRelaxedPrimalInfeasibilities_(0.0),
2446	acceptablePivot_(1.0e-8),
2447	lower_(NULL),
2448	rowLowerWork_(NULL),
2449	columnLowerWork_(NULL),
2450	upper_(NULL),
2451	rowUpperWork_(NULL),
2452	columnUpperWork_(NULL),
2453	cost_(NULL),
2454	rowObjectiveWork_(NULL),
2455	objectiveWork_(NULL),
2456	sequenceIn_(-1),
2457	directionIn_(-1),
2458	sequenceOut_(-1),
2459	directionOut_(-1),
2460	pivotRow_(-1),
2461	lastGoodIteration_(-100),
2462	dj_(NULL),
2463	rowReducedCost_(NULL),
2464	reducedCostWork_(NULL),
2465	solution_(NULL),
2466	rowActivityWork_(NULL),
2467	columnActivityWork_(NULL),
2468	numberDualInfeasibilities_(0),
2469	numberDualInfeasibilitiesWithoutFree_(0),
2470	numberPrimalInfeasibilities_(100),
2471	numberRefinements_(0),
2472	pivotVariable_(NULL),
2473	factorization_(NULL),
2474	savedSolution_(NULL),
2475	numberTimesOptimal_(0),
2476	disasterArea_(NULL),
2477	changeMade_(1),
2478	algorithm_(0),
2479	forceFactorization_(-1),
2480	perturbation_(100),
2481	nonLinearCost_(NULL),
2482	lastBadIteration_(-999999),
2483	lastFlaggedIteration_(-999999),
2484	numberFake_(0),
2485	numberChanged_(0),
2486	progressFlag_(0),
2487	firstFree_(-1),
2488	numberExtraRows_(0),
2489	maximumBasic_(0),
2490	dontFactorizePivots_(0),
2491	incomingInfeasibility_(1.0),
2492	allowedInfeasibility_(10.0),
2493	automaticScale_(0),
2494	maximumPerturbationSize_(0),
2495	perturbationArray_(NULL),
2496	baseModel_(NULL)
2497	#ifdef ABC_INHERIT
2498	,abcSimplex_(NULL),
2499	abcState_(0)
2500	#endif
2501	{
2502	int i;
2503	for (i = 0; i < 6; i++) {
2504	rowArray_[i] = NULL;
2505	columnArray_[i] = NULL;
2506	}
2507	for (i = 0; i < 4; i++) {
2508	spareIntArray_[i] = 0;
2509	spareDoubleArray_[i] = 0.0;
2510	}
2511	saveStatus_ = NULL;
2512	// get an empty factorization so we can set tolerances etc
2513	getEmptyFactorization();
2514	// say Steepest pricing
2515	dualRowPivot_ = new ClpDualRowSteepest();
2516	// say Steepest pricing
2517	primalColumnPivot_ = new ClpPrimalColumnSteepest();
2518	solveType_ = 1; // say simplex based life form
2519
2520	}
2521	// Assignment operator. This copies the data
2522	ClpSimplex &
2523	ClpSimplex::operator=(const ClpSimplex & rhs)
2524	{
2525	if (this != &rhs) {
2526	gutsOfDelete(0);
2527	delete nonLinearCost_;
2528	nonLinearCost_ = NULL;
2529	ClpModel::operator=(rhs);
2530	gutsOfCopy(rhs);
2531	}
2532	return *this;
2533	}
2534	void
2535	ClpSimplex::gutsOfCopy(const ClpSimplex & rhs)
2536	{
2537	assert (numberRows_ == rhs.numberRows_);
2538	assert (numberColumns_ == rhs.numberColumns_);
2539	numberExtraRows_ = rhs.numberExtraRows_;
2540	maximumBasic_ = rhs.maximumBasic_;
2541	dontFactorizePivots_ = rhs.dontFactorizePivots_;
2542	int numberRows2 = numberRows_ + numberExtraRows_;
2543	moreSpecialOptions_ = rhs.moreSpecialOptions_;
2544	if ((whatsChanged_ & 1) != 0) {
2545	int numberTotal = numberColumns_ + numberRows2;
2546	if ((specialOptions_ & 65536) != 0 && maximumRows_ >= 0) {
2547	assert (maximumInternalRows_ >= numberRows2);
2548	assert (maximumInternalColumns_ >= numberColumns_);
2549	numberTotal = 2 * (maximumInternalColumns_ + maximumInternalRows_);
2550	}
2551	lower_ = ClpCopyOfArray(rhs.lower_, numberTotal);
2552	rowLowerWork_ = lower_ + numberColumns_;
2553	columnLowerWork_ = lower_;
2554	upper_ = ClpCopyOfArray(rhs.upper_, numberTotal);
2555	rowUpperWork_ = upper_ + numberColumns_;
2556	columnUpperWork_ = upper_;
2557	cost_ = ClpCopyOfArray(rhs.cost_, numberTotal);
2558	objectiveWork_ = cost_;
2559	rowObjectiveWork_ = cost_ + numberColumns_;
2560	dj_ = ClpCopyOfArray(rhs.dj_, numberTotal);
2561	if (dj_) {
2562	reducedCostWork_ = dj_;
2563	rowReducedCost_ = dj_ + numberColumns_;
2564	}
2565	solution_ = ClpCopyOfArray(rhs.solution_, numberTotal);
2566	if (solution_) {
2567	columnActivityWork_ = solution_;
2568	rowActivityWork_ = solution_ + numberColumns_;
2569	}
2570	if (rhs.pivotVariable_) {
2571	pivotVariable_ = new int[numberRows2];
2572	CoinMemcpyN ( rhs.pivotVariable_, numberRows2 , pivotVariable_);
2573	} else {
2574	pivotVariable_ = NULL;
2575	}
2576	savedSolution_ = ClpCopyOfArray(rhs.savedSolution_, numberTotal);
2577	int i;
2578	for (i = 0; i < 6; i++) {
2579	rowArray_[i] = NULL;
2580	if (rhs.rowArray_[i])
2581	rowArray_[i] = new CoinIndexedVector(*rhs.rowArray_[i]);
2582	columnArray_[i] = NULL;
2583	if (rhs.columnArray_[i])
2584	columnArray_[i] = new CoinIndexedVector(*rhs.columnArray_[i]);
2585	}
2586	if (rhs.saveStatus_) {
2587	saveStatus_ = ClpCopyOfArray( rhs.saveStatus_, numberTotal);
2588	}
2589	} else {
2590	lower_ = NULL;
2591	rowLowerWork_ = NULL;
2592	columnLowerWork_ = NULL;
2593	upper_ = NULL;
2594	rowUpperWork_ = NULL;
2595	columnUpperWork_ = NULL;
2596	cost_ = NULL;
2597	objectiveWork_ = NULL;
2598	rowObjectiveWork_ = NULL;
2599	dj_ = NULL;
2600	reducedCostWork_ = NULL;
2601	rowReducedCost_ = NULL;
2602	solution_ = NULL;
2603	columnActivityWork_ = NULL;
2604	rowActivityWork_ = NULL;
2605	pivotVariable_ = NULL;
2606	savedSolution_ = NULL;
2607	int i;
2608	for (i = 0; i < 6; i++) {
2609	rowArray_[i] = NULL;
2610	columnArray_[i] = NULL;
2611	}
2612	saveStatus_ = NULL;
2613	}
2614	if (rhs.factorization_) {
2615	setFactorization(*rhs.factorization_);
2616	} else {
2617	delete factorization_;
2618	factorization_ = NULL;
2619	}
2620	bestPossibleImprovement_ = rhs.bestPossibleImprovement_;
2621	columnPrimalSequence_ = rhs.columnPrimalSequence_;
2622	zeroTolerance_ = rhs.zeroTolerance_;
2623	rowPrimalSequence_ = rhs.rowPrimalSequence_;
2624	bestObjectiveValue_ = rhs.bestObjectiveValue_;
2625	baseIteration_ = rhs.baseIteration_;
2626	primalToleranceToGetOptimal_ = rhs.primalToleranceToGetOptimal_;
2627	largeValue_ = rhs.largeValue_;
2628	largestPrimalError_ = rhs.largestPrimalError_;
2629	largestDualError_ = rhs.largestDualError_;
2630	alphaAccuracy_ = rhs.alphaAccuracy_;
2631	dualBound_ = rhs.dualBound_;
2632	alpha_ = rhs.alpha_;
2633	theta_ = rhs.theta_;
2634	lowerIn_ = rhs.lowerIn_;
2635	valueIn_ = rhs.valueIn_;
2636	upperIn_ = rhs.upperIn_;
2637	dualIn_ = rhs.dualIn_;
2638	sequenceIn_ = rhs.sequenceIn_;
2639	directionIn_ = rhs.directionIn_;
2640	lowerOut_ = rhs.lowerOut_;
2641	valueOut_ = rhs.valueOut_;
2642	upperOut_ = rhs.upperOut_;
2643	dualOut_ = rhs.dualOut_;
2644	sequenceOut_ = rhs.sequenceOut_;
2645	directionOut_ = rhs.directionOut_;
2646	pivotRow_ = rhs.pivotRow_;
2647	lastGoodIteration_ = rhs.lastGoodIteration_;
2648	numberRefinements_ = rhs.numberRefinements_;
2649	dualTolerance_ = rhs.dualTolerance_;
2650	primalTolerance_ = rhs.primalTolerance_;
2651	sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_;
2652	numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
2653	numberDualInfeasibilitiesWithoutFree_ =
2654	rhs.numberDualInfeasibilitiesWithoutFree_;
2655	sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
2656	numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
2657	dualRowPivot_ = rhs.dualRowPivot_->clone(true);
2658	dualRowPivot_->setModel(this);
2659	primalColumnPivot_ = rhs.primalColumnPivot_->clone(true);
2660	primalColumnPivot_->setModel(this);
2661	numberTimesOptimal_ = rhs.numberTimesOptimal_;
2662	disasterArea_ = NULL;
2663	changeMade_ = rhs.changeMade_;
2664	algorithm_ = rhs.algorithm_;
2665	forceFactorization_ = rhs.forceFactorization_;
2666	perturbation_ = rhs.perturbation_;
2667	infeasibilityCost_ = rhs.infeasibilityCost_;
2668	lastBadIteration_ = rhs.lastBadIteration_;
2669	lastFlaggedIteration_ = rhs.lastFlaggedIteration_;
2670	numberFake_ = rhs.numberFake_;
2671	numberChanged_ = rhs.numberChanged_;
2672	progressFlag_ = rhs.progressFlag_;
2673	firstFree_ = rhs.firstFree_;
2674	incomingInfeasibility_ = rhs.incomingInfeasibility_;
2675	allowedInfeasibility_ = rhs.allowedInfeasibility_;
2676	automaticScale_ = rhs.automaticScale_;
2677	#ifdef ABC_INHERIT
2678	abcSimplex_ = NULL;
2679	abcState_ = rhs.abcState_;
2680	#endif
2681	maximumPerturbationSize_ = rhs.maximumPerturbationSize_;
2682	if (maximumPerturbationSize_ && maximumPerturbationSize_ >= 2 * numberColumns_) {
2683	perturbationArray_ = CoinCopyOfArray(rhs.perturbationArray_,
2684	maximumPerturbationSize_);
2685	} else {
2686	maximumPerturbationSize_ = 0;
2687	perturbationArray_ = NULL;
2688	}
2689	if (rhs.baseModel_) {
2690	baseModel_ = new ClpSimplex(*rhs.baseModel_);
2691	} else {
2692	baseModel_ = NULL;
2693	}
2694	progress_ = rhs.progress_;
2695	for (int i = 0; i < 4; i++) {
2696	spareIntArray_[i] = rhs.spareIntArray_[i];
2697	spareDoubleArray_[i] = rhs.spareDoubleArray_[i];
2698	}
2699	sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
2700	sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
2701	acceptablePivot_ = rhs.acceptablePivot_;
2702	if (rhs.nonLinearCost_ != NULL)
2703	nonLinearCost_ = new ClpNonLinearCost(*rhs.nonLinearCost_);
2704	else
2705	nonLinearCost_ = NULL;
2706	solveType_ = rhs.solveType_;
2707	eventHandler_->setSimplex(this);
2708	}
2709	// type == 0 do everything, most + pivot data, 2 factorization data as well
2710	void
2711	ClpSimplex::gutsOfDelete(int type)
2712	{
2713	if (!type \|\| (specialOptions_ & 65536) == 0) {
2714	maximumInternalColumns_ = -1;
2715	maximumInternalRows_ = -1;
2716	delete [] lower_;
2717	lower_ = NULL;
2718	rowLowerWork_ = NULL;
2719	columnLowerWork_ = NULL;
2720	delete [] upper_;
2721	upper_ = NULL;
2722	rowUpperWork_ = NULL;
2723	columnUpperWork_ = NULL;
2724	delete [] cost_;
2725	cost_ = NULL;
2726	objectiveWork_ = NULL;
2727	rowObjectiveWork_ = NULL;
2728	delete [] dj_;
2729	dj_ = NULL;
2730	reducedCostWork_ = NULL;
2731	rowReducedCost_ = NULL;
2732	delete [] solution_;
2733	solution_ = NULL;
2734	rowActivityWork_ = NULL;
2735	columnActivityWork_ = NULL;
2736	delete [] savedSolution_;
2737	savedSolution_ = NULL;
2738	}
2739	if ((specialOptions_ & 2) == 0) {
2740	delete nonLinearCost_;
2741	nonLinearCost_ = NULL;
2742	}
2743	int i;
2744	if ((specialOptions_ & 65536) == 0) {
2745	for (i = 0; i < 6; i++) {
2746	delete rowArray_[i];
2747	rowArray_[i] = NULL;
2748	delete columnArray_[i];
2749	columnArray_[i] = NULL;
2750	}
2751	}
2752	delete [] saveStatus_;
2753	saveStatus_ = NULL;
2754	if (type != 1) {
2755	delete rowCopy_;
2756	rowCopy_ = NULL;
2757	}
2758	if (!type) {
2759	// delete everything
2760	setEmptyFactorization();
2761	delete [] pivotVariable_;
2762	pivotVariable_ = NULL;
2763	delete dualRowPivot_;
2764	dualRowPivot_ = NULL;
2765	delete primalColumnPivot_;
2766	primalColumnPivot_ = NULL;
2767	delete baseModel_;
2768	baseModel_ = NULL;
2769	delete [] perturbationArray_;
2770	perturbationArray_ = NULL;
2771	maximumPerturbationSize_ = 0;
2772	} else {
2773	// delete any size information in methods
2774	if (type > 1) {
2775	//assert (factorization_);
2776	if (factorization_)
2777	factorization_->clearArrays();
2778	delete [] pivotVariable_;
2779	pivotVariable_ = NULL;
2780	}
2781	dualRowPivot_->clearArrays();
2782	primalColumnPivot_->clearArrays();
2783	}
2784	}
2785	// This sets largest infeasibility and most infeasible
2786	void
2787	ClpSimplex::checkPrimalSolution(const double * rowActivities,
2788	const double * columnActivities)
2789	{
2790	double * solution;
2791	int iRow, iColumn;
2792
2793	objectiveValue_ = 0.0;
2794	// now look at primal solution
2795	solution = rowActivityWork_;
2796	sumPrimalInfeasibilities_ = 0.0;
2797	numberPrimalInfeasibilities_ = 0;
2798	double primalTolerance = primalTolerance_;
2799	double relaxedTolerance = primalTolerance_;
2800	// we can't really trust infeasibilities if there is primal error
2801	double error = CoinMin(1.0e-2, largestPrimalError_);
2802	// allow tolerance at least slightly bigger than standard
2803	relaxedTolerance = relaxedTolerance + error;
2804	sumOfRelaxedPrimalInfeasibilities_ = 0.0;
2805	for (iRow = 0; iRow < numberRows_; iRow++) {
2806	//assert (fabs(solution[iRow])<1.0e15\|\|getRowStatus(iRow) == basic);
2807	double infeasibility = 0.0;
2808	objectiveValue_ += solution[iRow] * rowObjectiveWork_[iRow];
2809	if (solution[iRow] > rowUpperWork_[iRow]) {
2810	infeasibility = solution[iRow] - rowUpperWork_[iRow];
2811	} else if (solution[iRow] < rowLowerWork_[iRow]) {
2812	infeasibility = rowLowerWork_[iRow] - solution[iRow];
2813	}
2814	if (infeasibility > primalTolerance) {
2815	sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2816	if (infeasibility > relaxedTolerance)
2817	sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2818	numberPrimalInfeasibilities_ ++;
2819	}
2820	infeasibility = fabs(rowActivities[iRow] - solution[iRow]);
2821	}
2822	// Check any infeasibilities from dynamic rows
2823	matrix_->primalExpanded(this, 2);
2824	solution = columnActivityWork_;
2825	if (!matrix_->rhsOffset(this)) {
2826	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2827	//assert (fabs(solution[iColumn])<1.0e15\|\|getColumnStatus(iColumn) == basic);
2828	double infeasibility = 0.0;
2829	objectiveValue_ += objectiveWork_[iColumn] * solution[iColumn];
2830	if (solution[iColumn] > columnUpperWork_[iColumn]) {
2831	infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
2832	} else if (solution[iColumn] < columnLowerWork_[iColumn]) {
2833	infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
2834	}
2835	if (infeasibility > primalTolerance) {
2836	sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2837	if (infeasibility > relaxedTolerance)
2838	sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2839	numberPrimalInfeasibilities_ ++;
2840	}
2841	infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
2842	}
2843	} else {
2844	// as we are using effective rhs we only check basics
2845	// But we do need to get objective
2846	objectiveValue_ += innerProduct(objectiveWork_, numberColumns_, solution);
2847	for (int j = 0; j < numberRows_; j++) {
2848	int iColumn = pivotVariable_[j];
2849	//assert (fabs(solution[iColumn])<1.0e15\|\|getColumnStatus(iColumn) == basic);
2850	double infeasibility = 0.0;
2851	if (solution[iColumn] > columnUpperWork_[iColumn]) {
2852	infeasibility = solution[iColumn] - columnUpperWork_[iColumn];
2853	} else if (solution[iColumn] < columnLowerWork_[iColumn]) {
2854	infeasibility = columnLowerWork_[iColumn] - solution[iColumn];
2855	}
2856	if (infeasibility > primalTolerance) {
2857	sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
2858	if (infeasibility > relaxedTolerance)
2859	sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedTolerance;
2860	numberPrimalInfeasibilities_ ++;
2861	}
2862	infeasibility = fabs(columnActivities[iColumn] - solution[iColumn]);
2863	}
2864	}
2865	objectiveValue_ += objective_->nonlinearOffset();
2866	objectiveValue_ /= (objectiveScale_ * rhsScale_);
2867	}
2868	void
2869	ClpSimplex::checkDualSolution()
2870	{
2871
2872	int iRow, iColumn;
2873	sumDualInfeasibilities_ = 0.0;
2874	numberDualInfeasibilities_ = 0;
2875	numberDualInfeasibilitiesWithoutFree_ = 0;
2876	if (matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) {
2877	// pretend we found dual infeasibilities
2878	sumOfRelaxedDualInfeasibilities_ = 1.0;
2879	sumDualInfeasibilities_ = 1.0;
2880	numberDualInfeasibilities_ = 1;
2881	return;
2882	}
2883	int firstFreePrimal = -1;
2884	int firstFreeDual = -1;
2885	int numberSuperBasicWithDj = 0;
2886	bestPossibleImprovement_ = 0.0;
2887	// we can't really trust infeasibilities if there is dual error
2888	double error = CoinMin(1.0e-2, largestDualError_);
2889	// allow tolerance at least slightly bigger than standard
2890	double relaxedTolerance = dualTolerance_ + error;
2891	// allow bigger tolerance for possible improvement
2892	double possTolerance = 5.0 * relaxedTolerance;
2893	sumOfRelaxedDualInfeasibilities_ = 0.0;
2894
2895	// Check any djs from dynamic rows
2896	matrix_->dualExpanded(this, NULL, NULL, 3);
2897	numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_;
2898	objectiveValue_ = 0.0;
2899	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2900	objectiveValue_ += objectiveWork_[iColumn] * columnActivityWork_[iColumn];
2901	if (getColumnStatus(iColumn) != basic && !flagged(iColumn)) {
2902	// not basic
2903	double distanceUp = columnUpperWork_[iColumn] -
2904	columnActivityWork_[iColumn];
2905	double distanceDown = columnActivityWork_[iColumn] -
2906	columnLowerWork_[iColumn];
2907	if (distanceUp > primalTolerance_) {
2908	double value = reducedCostWork_[iColumn];
2909	// Check if "free"
2910	if (distanceDown > primalTolerance_) {
2911	if (fabs(value) > 1.0e2 * relaxedTolerance) {
2912	numberSuperBasicWithDj++;
2913	if (firstFreeDual < 0)
2914	firstFreeDual = iColumn;
2915	}
2916	if (firstFreePrimal < 0)
2917	firstFreePrimal = iColumn;
2918	}
2919	// should not be negative
2920	if (value < 0.0) {
2921	value = - value;
2922	if (value > dualTolerance_) {
2923	if (getColumnStatus(iColumn) != isFree) {
2924	numberDualInfeasibilitiesWithoutFree_ ++;
2925	sumDualInfeasibilities_ += value - dualTolerance_;
2926	if (value > possTolerance)
2927	bestPossibleImprovement_ += CoinMin(distanceUp, 1.0e10) * value;
2928	if (value > relaxedTolerance)
2929	sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2930	numberDualInfeasibilities_ ++;
2931	} else {
2932	// free so relax a lot
2933	value *= 0.01;
2934	if (value > dualTolerance_) {
2935	sumDualInfeasibilities_ += value - dualTolerance_;
2936	if (value > possTolerance)
2937	bestPossibleImprovement_ = 1.0e100;
2938	if (value > relaxedTolerance)
2939	sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2940	numberDualInfeasibilities_ ++;
2941	}
2942	}
2943	}
2944	}
2945	}
2946	if (distanceDown > primalTolerance_) {
2947	double value = reducedCostWork_[iColumn];
2948	// should not be positive
2949	if (value > 0.0) {
2950	if (value > dualTolerance_) {
2951	sumDualInfeasibilities_ += value - dualTolerance_;
2952	if (value > possTolerance)
2953	bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
2954	if (value > relaxedTolerance)
2955	sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2956	numberDualInfeasibilities_ ++;
2957	if (getColumnStatus(iColumn) != isFree)
2958	numberDualInfeasibilitiesWithoutFree_ ++;
2959	// maybe we can make feasible by increasing tolerance
2960	}
2961	}
2962	}
2963	}
2964	}
2965	for (iRow = 0; iRow < numberRows_; iRow++) {
2966	objectiveValue_ += rowActivityWork_[iRow] * rowObjectiveWork_[iRow];
2967	if (getRowStatus(iRow) != basic && !flagged(iRow + numberColumns_)) {
2968	// not basic
2969	double distanceUp = rowUpperWork_[iRow] - rowActivityWork_[iRow];
2970	double distanceDown = rowActivityWork_[iRow] - rowLowerWork_[iRow];
2971	if (distanceUp > primalTolerance_) {
2972	double value = rowReducedCost_[iRow];
2973	// Check if "free"
2974	if (distanceDown > primalTolerance_) {
2975	if (fabs(value) > 1.0e2 * relaxedTolerance) {
2976	numberSuperBasicWithDj++;
2977	if (firstFreeDual < 0)
2978	firstFreeDual = iRow + numberColumns_;
2979	}
2980	if (firstFreePrimal < 0)
2981	firstFreePrimal = iRow + numberColumns_;
2982	}
2983	// should not be negative
2984	if (value < 0.0) {
2985	value = - value;
2986	if (value > dualTolerance_) {
2987	sumDualInfeasibilities_ += value - dualTolerance_;
2988	if (value > possTolerance)
2989	bestPossibleImprovement_ += value * CoinMin(distanceUp, 1.0e10);
2990	if (value > relaxedTolerance)
2991	sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
2992	numberDualInfeasibilities_ ++;
2993	if (getRowStatus(iRow) != isFree)
2994	numberDualInfeasibilitiesWithoutFree_ ++;
2995	}
2996	}
2997	}
2998	if (distanceDown > primalTolerance_) {
2999	double value = rowReducedCost_[iRow];
3000	// should not be positive
3001	if (value > 0.0) {
3002	if (value > dualTolerance_) {
3003	sumDualInfeasibilities_ += value - dualTolerance_;
3004	if (value > possTolerance)
3005	bestPossibleImprovement_ += value * CoinMin(distanceDown, 1.0e10);
3006	if (value > relaxedTolerance)
3007	sumOfRelaxedDualInfeasibilities_ += value - relaxedTolerance;
3008	numberDualInfeasibilities_ ++;
3009	if (getRowStatus(iRow) != isFree)
3010	numberDualInfeasibilitiesWithoutFree_ ++;
3011	// maybe we can make feasible by increasing tolerance
3012	}
3013	}
3014	}
3015	}
3016	}
3017	if (algorithm_ < 0 && firstFreeDual >= 0) {
3018	// dual
3019	firstFree_ = firstFreeDual;
3020	} else if (numberSuperBasicWithDj \|\|
3021	(progress_.lastIterationNumber(0) <= 0)) {
3022	firstFree_ = firstFreePrimal;
3023	}
3024	objectiveValue_ += objective_->nonlinearOffset();
3025	objectiveValue_ /= (objectiveScale_ * rhsScale_);
3026	}
3027	/* This sets sum and number of infeasibilities (Dual and Primal) */
3028	void
3029	ClpSimplex::checkBothSolutions()
3030	{
3031	if ((matrix_->skipDualCheck() && algorithm_ > 0 && problemStatus_ == -2) \|\|
3032	matrix_->rhsOffset(this)) {
3033	// Say may be free or superbasic
3034	moreSpecialOptions_ &= ~8;
3035	// old way
3036	checkPrimalSolution(rowActivityWork_, columnActivityWork_);
3037	checkDualSolution();
3038	return;
3039	}
3040	int iSequence;
3041	assert (dualTolerance_ > 0.0 && dualTolerance_ < 1.0e10);
3042	assert (primalTolerance_ > 0.0 && primalTolerance_ < 1.0e10);
3043	objectiveValue_ = 0.0;
3044	sumPrimalInfeasibilities_ = 0.0;
3045	numberPrimalInfeasibilities_ = 0;
3046	double primalTolerance = primalTolerance_;
3047	double relaxedToleranceP = primalTolerance_;
3048	// we can't really trust infeasibilities if there is primal error
3049	double error = CoinMin(1.0e-2, CoinMax(largestPrimalError_,0.0*primalTolerance_));
3050	// allow tolerance at least slightly bigger than standard
3051	relaxedToleranceP = relaxedToleranceP + error;
3052	sumOfRelaxedPrimalInfeasibilities_ = 0.0;
3053	sumDualInfeasibilities_ = 0.0;
3054	numberDualInfeasibilities_ = 0;
3055	double dualTolerance = dualTolerance_;
3056	double relaxedToleranceD = dualTolerance;
3057	// we can't really trust infeasibilities if there is dual error
3058	error = CoinMin(1.0e-2, CoinMax(largestDualError_,5.0*dualTolerance_));
3059	// allow tolerance at least slightly bigger than standard
3060	relaxedToleranceD = relaxedToleranceD + error;
3061	// allow bigger tolerance for possible improvement
3062	double possTolerance = 5.0 * relaxedToleranceD;
3063	sumOfRelaxedDualInfeasibilities_ = 0.0;
3064	bestPossibleImprovement_ = 0.0;
3065
3066	// Check any infeasibilities from dynamic rows
3067	matrix_->primalExpanded(this, 2);
3068	// Check any djs from dynamic rows
3069	matrix_->dualExpanded(this, NULL, NULL, 3);
3070	int numberDualInfeasibilitiesFree = 0;
3071	int firstFreePrimal = -1;
3072	int firstFreeDual = -1;
3073	int numberSuperBasicWithDj = 0;
3074
3075	int numberTotal = numberRows_ + numberColumns_;
3076	// Say no free or superbasic
3077	moreSpecialOptions_ \|= 8;
3078	//#define PRINT_INFEAS
3079	#ifdef PRINT_INFEAS
3080	int seqInf[10];
3081	#endif
3082	for (iSequence = 0; iSequence < numberTotal; iSequence++) {
3083	double value = solution_[iSequence];
3084	#ifdef COIN_DEBUG
3085	if (fabs(value) > 1.0e20)
3086	printf("%d values %g %g %g - status %d\n", iSequence, lower_[iSequence],
3087	solution_[iSequence], upper_[iSequence], status_[iSequence]);
3088	#endif
3089	objectiveValue_ += value * cost_[iSequence];
3090	double distanceUp = upper_[iSequence] - value;
3091	double distanceDown = value - lower_[iSequence];
3092	if (distanceUp < -primalTolerance) {
3093	double infeasibility = -distanceUp;
3094	sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
3095	if (infeasibility > relaxedToleranceP)
3096	sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
3097	#ifdef PRINT_INFEAS
3098	if (numberPrimalInfeasibilities_<10) {
3099	seqInf[numberPrimalInfeasibilities_]=iSequence;
3100	}
3101	#endif
3102	numberPrimalInfeasibilities_ ++;
3103	} else if (distanceDown < -primalTolerance) {
3104	double infeasibility = -distanceDown;
3105	sumPrimalInfeasibilities_ += infeasibility - primalTolerance_;
3106	if (infeasibility > relaxedToleranceP)
3107	sumOfRelaxedPrimalInfeasibilities_ += infeasibility - relaxedToleranceP;
3108	#ifdef PRINT_INFEAS
3109	if (numberPrimalInfeasibilities_<10) {
3110	seqInf[numberPrimalInfeasibilities_]=iSequence;
3111	}
3112	#endif
3113	numberPrimalInfeasibilities_ ++;
3114	} else {
3115	// feasible (so could be free)
3116	if (getStatus(iSequence) != basic && !flagged(iSequence)) {
3117	// not basic
3118	double djValue = dj_[iSequence];
3119	if (distanceDown < primalTolerance) {
3120	if (distanceUp > primalTolerance && djValue < -dualTolerance) {
3121	sumDualInfeasibilities_ -= djValue + dualTolerance;
3122	if (djValue < -possTolerance)
3123	bestPossibleImprovement_ -= distanceUp * djValue;
3124	if (djValue < -relaxedToleranceD)
3125	sumOfRelaxedDualInfeasibilities_ -= djValue + relaxedToleranceD;
3126	numberDualInfeasibilities_ ++;
3127	}
3128	} else if (distanceUp < primalTolerance) {
3129	if (djValue > dualTolerance) {
3130	sumDualInfeasibilities_ += djValue - dualTolerance;
3131	if (djValue > possTolerance)
3132	bestPossibleImprovement_ += distanceDown * djValue;
3133	if (djValue > relaxedToleranceD)
3134	sumOfRelaxedDualInfeasibilities_ += djValue - relaxedToleranceD;
3135	numberDualInfeasibilities_ ++;
3136	}
3137	} else {
3138	// may be free
3139	// Say free or superbasic
3140	moreSpecialOptions_ &= ~8;
3141	djValue *= 0.01;
3142	if (fabs(djValue) > dualTolerance) {
3143	if (getStatus(iSequence) == isFree)
3144	numberDualInfeasibilitiesFree++;
3145	sumDualInfeasibilities_ += fabs(djValue) - dualTolerance;
3146	bestPossibleImprovement_ = 1.0e100;
3147	numberDualInfeasibilities_ ++;
3148	if (fabs(djValue) > relaxedToleranceD) {
3149	sumOfRelaxedDualInfeasibilities_ += value - relaxedToleranceD;
3150	numberSuperBasicWithDj++;
3151	if (firstFreeDual < 0)
3152	firstFreeDual = iSequence;
3153	}
3154	}
3155	if (firstFreePrimal < 0)
3156	firstFreePrimal = iSequence;
3157	}
3158	}
3159	}
3160	}
3161	objectiveValue_ += objective_->nonlinearOffset();
3162	objectiveValue_ /= (objectiveScale_ * rhsScale_);
3163	numberDualInfeasibilitiesWithoutFree_ = numberDualInfeasibilities_ -
3164	numberDualInfeasibilitiesFree;
3165	#ifdef PRINT_INFEAS
3166	if (numberPrimalInfeasibilities_<=10) {
3167	printf("---------------start-----------\n");
3168	if (!rowScale_) {
3169	for (int i=0;i<numberPrimalInfeasibilities_;i++) {
3170	int iSeq = seqInf[i];
3171	double infeas;
3172	if (solution_[iSeq]<lower_[iSeq])
3173	infeas = lower_[iSeq]-solution_[iSeq];
3174	else
3175	infeas = solution_[iSeq]-upper_[iSeq];
3176	if (iSeq<numberColumns_) {
3177	printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g\n",
3178	iSeq,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas);
3179	} else {
3180	printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g\n",
3181	iSeq-numberColumns_,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas);
3182	}
3183	}
3184	} else {
3185	for (int i=0;i<numberPrimalInfeasibilities_;i++) {
3186	int iSeq = seqInf[i];
3187	double infeas;
3188	if (solution_[iSeq]<lower_[iSeq])
3189	infeas = lower_[iSeq]-solution_[iSeq];
3190	else
3191	infeas = solution_[iSeq]-upper_[iSeq];
3192	double unscaled = infeas;
3193	if (iSeq<numberColumns_) {
3194	unscaled *= columnScale_[iSeq];
3195	printf("INF C%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
3196	iSeq,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas,unscaled);
3197	} else {
3198	unscaled /= rowScale_[iSeq-numberColumns_];
3199	printf("INF R%d %.10g <= %.10g <= %.10g - infeas %g - unscaled %g\n",
3200	iSeq-numberColumns_,lower_[iSeq],solution_[iSeq],upper_[iSeq],infeas,unscaled);
3201	}
3202	}
3203	}
3204	}
3205	#endif
3206	if (algorithm_ < 0 && firstFreeDual >= 0) {
3207	// dual
3208	firstFree_ = firstFreeDual;
3209	} else if (numberSuperBasicWithDj \|\|
3210	(progress_.lastIterationNumber(0) <= 0)) {
3211	firstFree_ = firstFreePrimal;
3212	}
3213	}
3214	/* Adds multiple of a column into an array */
3215	void
3216	ClpSimplex::add(double * array,
3217	int sequence, double multiplier) const
3218	{
3219	if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3220	//slack
3221	array [sequence-numberColumns_] -= multiplier;
3222	} else {
3223	// column
3224	matrix_->add(this, array, sequence, multiplier);
3225	}
3226	}
3227	/*
3228	Unpacks one column of the matrix into indexed array
3229	*/
3230	void
3231	ClpSimplex::unpack(CoinIndexedVector * rowArray) const
3232	{
3233	rowArray->clear();
3234	if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
3235	//slack
3236	rowArray->insert(sequenceIn_ - numberColumns_, -1.0);
3237	} else {
3238	// column
3239	matrix_->unpack(this, rowArray, sequenceIn_);
3240	}
3241	}
3242	void
3243	ClpSimplex::unpack(CoinIndexedVector * rowArray, int sequence) const
3244	{
3245	rowArray->clear();
3246	if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3247	//slack
3248	rowArray->insert(sequence - numberColumns_, -1.0);
3249	} else {
3250	// column
3251	matrix_->unpack(this, rowArray, sequence);
3252	}
3253	}
3254	/*
3255	Unpacks one column of the matrix into indexed array
3256	*/
3257	void
3258	ClpSimplex::unpackPacked(CoinIndexedVector * rowArray)
3259	{
3260	rowArray->clear();
3261	if (sequenceIn_ >= numberColumns_ && sequenceIn_ < numberColumns_ + numberRows_) {
3262	//slack
3263	int * index = rowArray->getIndices();
3264	double * array = rowArray->denseVector();
3265	array[0] = -1.0;
3266	index[0] = sequenceIn_ - numberColumns_;
3267	rowArray->setNumElements(1);
3268	rowArray->setPackedMode(true);
3269	} else {
3270	// column
3271	matrix_->unpackPacked(this, rowArray, sequenceIn_);
3272	}
3273	}
3274	void
3275	ClpSimplex::unpackPacked(CoinIndexedVector * rowArray, int sequence)
3276	{
3277	rowArray->clear();
3278	if (sequence >= numberColumns_ && sequence < numberColumns_ + numberRows_) {
3279	//slack
3280	int * index = rowArray->getIndices();
3281	double * array = rowArray->denseVector();
3282	array[0] = -1.0;
3283	index[0] = sequence - numberColumns_;
3284	rowArray->setNumElements(1);
3285	rowArray->setPackedMode(true);
3286	} else {
3287	// column
3288	matrix_->unpackPacked(this, rowArray, sequence);
3289	}
3290	}
3291	//static int x_gaps[4]={0,0,0,0};
3292	//static int scale_times[]={0,0,0,0};
3293	bool
3294	ClpSimplex::createRim(int what, bool makeRowCopy, int startFinishOptions)
3295	{
3296	bool goodMatrix = true;
3297	int saveLevel = handler_->logLevel();
3298	spareIntArray_[0] = 0;
3299	if (!matrix_->canGetRowCopy())
3300	makeRowCopy = false; // switch off row copy if can't produce
3301	// Arrays will be there and correct size unless what is 63
3302	bool newArrays = (what == 63);
3303	// We may be restarting with same size
3304	bool keepPivots = false;
3305	if (startFinishOptions == -1) {
3306	startFinishOptions = 0;
3307	keepPivots = true;
3308	}
3309	bool oldMatrix = ((startFinishOptions & 4) != 0 && (whatsChanged_ & 1) != 0);
3310	if (what == 63) {
3311	pivotRow_ = -1;
3312	if (!status_)
3313	createStatus();
3314	if (oldMatrix)
3315	newArrays = false;
3316	if (problemStatus_ == 10) {
3317	handler_->setLogLevel(0); // switch off messages
3318	if (rowArray_[0]) {
3319	// stuff is still there
3320	oldMatrix = true;
3321	newArrays = false;
3322	keepPivots = true;
3323	for (int iRow = 0; iRow < 4; iRow++) {
3324	rowArray_[iRow]->clear();
3325	}
3326	for (int iColumn = 0; iColumn < 2; iColumn++) {
3327	columnArray_[iColumn]->clear();
3328	}
3329	}
3330	} else if (factorization_) {
3331	// match up factorization messages
3332	if (handler_->logLevel() < 3)
3333	factorization_->messageLevel(0);
3334	else
3335	factorization_->messageLevel(CoinMax(3, factorization_->messageLevel()));
3336	/* Faster to keep pivots rather than re-scan matrix. Matrix may have changed
3337	i.e. oldMatrix false but okay as long as same number rows and status array exists
3338	*/
3339	if ((startFinishOptions & 2) != 0 && factorization_->numberRows() == numberRows_ && status_)
3340	keepPivots = true;
3341	}
3342	numberExtraRows_ = matrix_->generalExpanded(this, 2, maximumBasic_);
3343	if (numberExtraRows_ && newArrays) {
3344	// make sure status array large enough
3345	assert (status_);
3346	int numberOld = numberRows_ + numberColumns_;
3347	int numberNew = numberRows_ + numberColumns_ + numberExtraRows_;
3348	unsigned char * newStatus = new unsigned char [numberNew];
3349	memset(newStatus + numberOld, 0, numberExtraRows_);
3350	CoinMemcpyN(status_, numberOld, newStatus);
3351	delete [] status_;
3352	status_ = newStatus;
3353	}
3354	}
3355	int numberRows2 = numberRows_ + numberExtraRows_;
3356	int numberTotal = numberRows2 + numberColumns_;
3357	if ((specialOptions_ & 65536) != 0) {
3358	assert (!numberExtraRows_);
3359	if (!cost_ \|\| numberRows2 > maximumInternalRows_ \|\|
3360	numberColumns_ > maximumInternalColumns_) {
3361	newArrays = true;
3362	keepPivots = false;
3363	COIN_DETAIL_PRINT(printf("createrim a %d rows, %d maximum rows %d maxinternal\n",
3364	numberRows_, maximumRows_, maximumInternalRows_));
3365	int oldMaximumRows = maximumInternalRows_;
3366	int oldMaximumColumns = maximumInternalColumns_;
3367	if (cost_) {
3368	if (numberRows2 > maximumInternalRows_)
3369	maximumInternalRows_ = numberRows2;
3370	if (numberColumns_ > maximumInternalColumns_)
3371	maximumInternalColumns_ = numberColumns_;
3372	} else {
3373	maximumInternalRows_ = numberRows2;
3374	maximumInternalColumns_ = numberColumns_;
3375	}
3376	//maximumRows_=CoinMax(maximumInternalRows_,maximumRows_);
3377	//maximumColumns_=CoinMax(maximumInternalColumns_,maximumColumns_);
3378	assert(maximumInternalRows_ == maximumRows_);
3379	assert(maximumInternalColumns_ == maximumColumns_);
3380	COIN_DETAIL_PRINT(printf("createrim b %d rows, %d maximum rows, %d maxinternal\n",
3381	numberRows_, maximumRows_, maximumInternalRows_));
3382	int numberTotal2 = (maximumInternalRows_ + maximumInternalColumns_) * 2;
3383	delete [] cost_;
3384	cost_ = new double[numberTotal2];
3385	delete [] lower_;
3386	delete [] upper_;
3387	lower_ = new double[numberTotal2];
3388	upper_ = new double[numberTotal2];
3389	delete [] dj_;
3390	dj_ = new double[numberTotal2];
3391	delete [] solution_;
3392	solution_ = new double[numberTotal2];
3393	// ***** should be non NULL but seems to be too much
3394	//printf("resize %d savedRowScale %x\n",maximumRows_,savedRowScale_);
3395	if (savedRowScale_) {
3396	assert (oldMaximumRows > 0);
3397	double * temp;
3398	temp = new double [4*maximumRows_];
3399	CoinFillN(temp, 4 * maximumRows_, 1.0);
3400	CoinMemcpyN(savedRowScale_, numberRows_, temp);
3401	CoinMemcpyN(savedRowScale_ + oldMaximumRows, numberRows_, temp + maximumRows_);
3402	CoinMemcpyN(savedRowScale_ + 2 * oldMaximumRows, numberRows_, temp + 2 * maximumRows_);
3403	CoinMemcpyN(savedRowScale_ + 3 * oldMaximumRows, numberRows_, temp + 3 * maximumRows_);
3404	delete [] savedRowScale_;
3405	savedRowScale_ = temp;
3406	temp = new double [4*maximumColumns_];
3407	CoinFillN(temp, 4 * maximumColumns_, 1.0);
3408	CoinMemcpyN(savedColumnScale_, numberColumns_, temp);
3409	CoinMemcpyN(savedColumnScale_ + oldMaximumColumns, numberColumns_, temp + maximumColumns_);
3410	CoinMemcpyN(savedColumnScale_ + 2 * oldMaximumColumns, numberColumns_, temp + 2 * maximumColumns_);
3411	CoinMemcpyN(savedColumnScale_ + 3 * oldMaximumColumns, numberColumns_, temp + 3 * maximumColumns_);
3412	delete [] savedColumnScale_;
3413	savedColumnScale_ = temp;
3414	}
3415	}
3416	}
3417	int i;
3418	bool doSanityCheck = true;
3419	if (what == 63) {
3420	// We may want to switch stuff off for speed
3421	if ((specialOptions_ & 256) != 0)
3422	makeRowCopy = false; // no row copy
3423	if ((specialOptions_ & 128) != 0)
3424	doSanityCheck = false; // no sanity check
3425	//check matrix
3426	if (!matrix_)
3427	matrix_ = new ClpPackedMatrix();
3428	int checkType = (doSanityCheck) ? 15 : 14;
3429	if (oldMatrix)
3430	checkType = 14;
3431	bool inCbcOrOther = (specialOptions_ & 0x03000000) != 0;
3432	if (inCbcOrOther)
3433	checkType -= 4; // don't check for duplicates
3434	if (!matrix_->allElementsInRange(this, smallElement_, 1.0e20, checkType)) {
3435	problemStatus_ = 4;
3436	secondaryStatus_ = 8;
3437	//goodMatrix= false;
3438	return false;
3439	}
3440	bool rowCopyIsScaled;
3441	if (makeRowCopy) {
3442	if(!oldMatrix \|\| !rowCopy_) {
3443	delete rowCopy_;
3444	// may return NULL if can't give row copy
3445	rowCopy_ = matrix_->reverseOrderedCopy();
3446	rowCopyIsScaled = false;
3447	} else {
3448	rowCopyIsScaled = true;
3449	}
3450	}
3451	#if 0
3452	if (what == 63) {
3453	int k = rowScale_ ? 1 : 0;
3454	if (oldMatrix)
3455	k += 2;
3456	scale_times[k]++;
3457	if ((scale_times[0] + scale_times[1] + scale_times[2] + scale_times[3]) % 1000 == 0)
3458	printf("scale counts %d %d %d %d\n",
3459	scale_times[0], scale_times[1], scale_times[2], scale_times[3]);
3460	}
3461	#endif
3462	// do scaling if needed
3463	if (!oldMatrix && scalingFlag_ < 0) {
3464	if (scalingFlag_ < 0 && rowScale_) {
3465	//if (handler_->logLevel()>0)
3466	printf("How did we get scalingFlag_ %d and non NULL rowScale_? - switching off scaling\n",
3467	scalingFlag_);
3468	scalingFlag_ = 0;
3469	}
3470	delete [] rowScale_;
3471	delete [] columnScale_;
3472	rowScale_ = NULL;
3473	columnScale_ = NULL;
3474	}
3475	inverseRowScale_ = NULL;
3476	inverseColumnScale_ = NULL;
3477	if (scalingFlag_ > 0 &&(specialOptions_ & 65536) != 0&&
3478	rowScale_&&rowScale_==savedRowScale_)
3479	rowScale_=NULL;
3480	if (scalingFlag_ > 0 && !rowScale_) {
3481	if ((specialOptions_ & 65536) != 0) {
3482	assert (!rowScale_);
3483	rowScale_ = savedRowScale_;
3484	columnScale_ = savedColumnScale_;
3485	// put back original
3486	if (savedRowScale_) {
3487	inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
3488	inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
3489	CoinMemcpyN(savedRowScale_ + 2 * maximumInternalRows_,
3490	numberRows2, savedRowScale_);
3491	CoinMemcpyN(savedRowScale_ + 3 * maximumInternalRows_,
3492	numberRows2, inverseRowScale_);
3493	CoinMemcpyN(savedColumnScale_ + 2 * maximumColumns_,
3494	numberColumns_, savedColumnScale_);
3495	CoinMemcpyN(savedColumnScale_ + 3 * maximumColumns_,
3496	numberColumns_, inverseColumnScale_);
3497	}
3498	}
3499	if (matrix_->scale(this))
3500	scalingFlag_ = -scalingFlag_; // not scaled after all
3501	if (rowScale_ && automaticScale_) {
3502	if (!savedRowScale_) {
3503	inverseRowScale_ = rowScale_ + numberRows2;
3504	inverseColumnScale_ = columnScale_ + numberColumns_;
3505	}
3506	// try automatic scaling
3507	double smallestObj = 1.0e100;
3508	double largestObj = 0.0;
3509	double largestRhs = 0.0;
3510	const double * obj = objective();
3511	for (i = 0; i < numberColumns_; i++) {
3512	double value = fabs(obj[i]);
3513	value *= columnScale_[i];
3514	if (value && columnLower_[i] != columnUpper_[i]) {
3515	smallestObj = CoinMin(smallestObj, value);
3516	largestObj = CoinMax(largestObj, value);
3517	}
3518	if (columnLower_[i] > 0.0 \|\| columnUpper_[i] < 0.0) {
3519	double scale = 1.0 * inverseColumnScale_[i];
3520	//printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3521	if (columnLower_[i] > 0)
3522	largestRhs = CoinMax(largestRhs, columnLower_[i] * scale);
3523	if (columnUpper_[i] < 0.0)
3524	largestRhs = CoinMax(largestRhs, -columnUpper_[i] * scale);
3525	}
3526	}
3527	for (i = 0; i < numberRows_; i++) {
3528	if (rowLower_[i] > 0.0 \|\| rowUpper_[i] < 0.0) {
3529	double scale = rowScale_[i];
3530	//printf("%d %g %g %g %g\n",i,scale,lower_[i],upper_[i],largestRhs);
3531	if (rowLower_[i] > 0)
3532	largestRhs = CoinMax(largestRhs, rowLower_[i] * scale);
3533	if (rowUpper_[i] < 0.0)
3534	largestRhs = CoinMax(largestRhs, -rowUpper_[i] * scale);
3535	}
3536	}
3537	COIN_DETAIL_PRINT(printf("small obj %g, large %g - rhs %g\n", smallestObj, largestObj, largestRhs));
3538	bool scalingDone = false;
3539	// look at element range
3540	double smallestNegative;
3541	double largestNegative;
3542	double smallestPositive;
3543	double largestPositive;
3544	matrix_->rangeOfElements(smallestNegative, largestNegative,
3545	smallestPositive, largestPositive);
3546	smallestPositive = CoinMin(fabs(smallestNegative), smallestPositive);
3547	largestPositive = CoinMax(fabs(largestNegative), largestPositive);
3548	if (largestObj) {
3549	double ratio = largestObj / smallestObj;
3550	double scale = 1.0;
3551	if (ratio < 1.0e8) {
3552	// reasonable
3553	if (smallestObj < 1.0e-4) {
3554	// may as well scale up
3555	scalingDone = true;
3556	scale = 1.0e-3 / smallestObj;
3557	} else if (largestObj < 1.0e6 \|\| (algorithm_ > 0 && largestObj < 1.0e-4 * infeasibilityCost_)) {
3558	//done=true;
3559	} else {
3560	scalingDone = true;
3561	if (algorithm_ < 0) {
3562	scale = 1.0e6 / largestObj;
3563	} else {
3564	scale = CoinMax(1.0e6, 1.0e-4 * infeasibilityCost_) / largestObj;
3565	}
3566	}
3567	} else if (ratio < 1.0e12) {
3568	// not so good
3569	if (smallestObj < 1.0e-7) {
3570	// may as well scale up
3571	scalingDone = true;
3572	scale = 1.0e-6 / smallestObj;
3573	} else if (largestObj < 1.0e7 \|\| (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3574	//done=true;
3575	} else {
3576	scalingDone = true;
3577	if (algorithm_ < 0) {
3578	scale = 1.0e7 / largestObj;
3579	} else {
3580	scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3581	}
3582	}
3583	} else {
3584	// Really nasty problem
3585	if (smallestObj < 1.0e-8) {
3586	// may as well scale up
3587	scalingDone = true;
3588	scale = 1.0e-7 / smallestObj;
3589	largestObj *= scale;
3590	}
3591	if (largestObj < 1.0e7 \|\| (algorithm_ > 0 && largestObj < 1.0e-3 * infeasibilityCost_)) {
3592	//done=true;
3593	} else {
3594	scalingDone = true;
3595	if (algorithm_ < 0) {
3596	scale = 1.0e7 / largestObj;
3597	} else {
3598	scale = CoinMax(1.0e7, 1.0e-3 * infeasibilityCost_) / largestObj;
3599	}
3600	}
3601	}
3602	objectiveScale_ = scale;
3603	}
3604	if (largestRhs > 1.0e12) {
3605	scalingDone = true;
3606	rhsScale_ = 1.0e9 / largestRhs;
3607	} else if (largestPositive > 1.0e-14 * smallestPositive && largestRhs > 1.0e6) {
3608	scalingDone = true;
3609	rhsScale_ = 1.0e6 / largestRhs;
3610	} else {
3611	rhsScale_ = 1.0;
3612	}
3613	if (scalingDone) {
3614	handler_->message(CLP_RIM_SCALE, messages_)
3615	<< objectiveScale_ << rhsScale_
3616	<< CoinMessageEol;
3617	}
3618	}
3619	} else if (makeRowCopy && scalingFlag_ > 0 && !rowCopyIsScaled) {
3620	matrix_->scaleRowCopy(this);
3621	}
3622	if (rowScale_ && !savedRowScale_) {
3623	inverseRowScale_ = rowScale_ + numberRows2;
3624	inverseColumnScale_ = columnScale_ + numberColumns_;
3625	}
3626	// See if we can try for faster row copy
3627	if (makeRowCopy && !oldMatrix) {
3628	ClpPackedMatrix* clpMatrix =
3629	dynamic_cast< ClpPackedMatrix*>(matrix_);
3630	if (clpMatrix && numberThreads_)
3631	clpMatrix->specialRowCopy(this, rowCopy_);
3632	if (clpMatrix)
3633	clpMatrix->specialColumnCopy(this);
3634	}
3635	}
3636	if (what == 63) {
3637	#if 0
3638	{
3639	x_gaps[0]++;
3640	ClpPackedMatrix* clpMatrix =
3641	dynamic_cast< ClpPackedMatrix*>(matrix_);
3642	if (clpMatrix) {
3643	if (!clpMatrix->getPackedMatrix()->hasGaps())
3644	x_gaps[1]++;
3645	if ((clpMatrix->flags() & 2) == 0)
3646	x_gaps[3]++;
3647	} else {
3648	x_gaps[2]++;
3649	}
3650	if ((x_gaps[0] % 1000) == 0)
3651	printf("create %d times, no gaps %d times - not clp %d times - flagged %d\n",
3652	x_gaps[0], x_gaps[1], x_gaps[2], x_gaps[3]);
3653	}
3654	#endif
3655	if (newArrays && (specialOptions_ & 65536) == 0) {
3656	delete [] cost_;
3657	cost_ = new double[2*numberTotal];
3658	delete [] lower_;
3659	delete [] upper_;
3660	lower_ = new double[numberTotal];
3661	upper_ = new double[numberTotal];
3662	delete [] dj_;
3663	dj_ = new double[numberTotal];
3664	delete [] solution_;
3665	solution_ = new double[numberTotal];
3666	}
3667	reducedCostWork_ = dj_;
3668	rowReducedCost_ = dj_ + numberColumns_;
3669	columnActivityWork_ = solution_;
3670	rowActivityWork_ = solution_ + numberColumns_;
3671	objectiveWork_ = cost_;
3672	rowObjectiveWork_ = cost_ + numberColumns_;
3673	rowLowerWork_ = lower_ + numberColumns_;
3674	columnLowerWork_ = lower_;
3675	rowUpperWork_ = upper_ + numberColumns_;
3676	columnUpperWork_ = upper_;
3677	}
3678	if ((what & 4) != 0) {
3679	double direction = optimizationDirection_ * objectiveScale_;
3680	const double * obj = objective();
3681	const double * rowScale = rowScale_;
3682	const double * columnScale = columnScale_;
3683	// and also scale by scale factors
3684	if (rowScale) {
3685	if (rowObjective_) {
3686	for (i = 0; i < numberRows_; i++)
3687	rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
3688	} else {
3689	memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3690	}
3691	// If scaled then do all columns later in one loop
3692	if (what != 63) {
3693	for (i = 0; i < numberColumns_; i++) {
3694	CoinAssert(fabs(obj[i]) < 1.0e25);
3695	objectiveWork_[i] = obj[i] * direction * columnScale[i];
3696	}
3697	}
3698	} else {
3699	if (rowObjective_) {
3700	for (i = 0; i < numberRows_; i++)
3701	rowObjectiveWork_[i] = rowObjective_[i] * direction;
3702	} else {
3703	memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
3704	}
3705	for (i = 0; i < numberColumns_; i++) {
3706	CoinAssert(fabs(obj[i]) < 1.0e25);
3707	objectiveWork_[i] = obj[i] * direction;
3708	}
3709	}
3710	}
3711	if ((what & 1) != 0) {
3712	const double * rowScale = rowScale_;
3713	// clean up any mismatches on infinity
3714	// and fix any variables with tiny gaps
3715	double primalTolerance = dblParam_[ClpPrimalTolerance];
3716	if(rowScale) {
3717	// If scaled then do all columns later in one loop
3718	if (what != 63) {
3719	const double * inverseScale = inverseColumnScale_;
3720	for (i = 0; i < numberColumns_; i++) {
3721	double multiplier = rhsScale_ * inverseScale[i];
3722	double lowerValue = columnLower_[i];
3723	double upperValue = columnUpper_[i];
3724	if (lowerValue > -1.0e20) {
3725	columnLowerWork_[i] = lowerValue * multiplier;
3726	if (upperValue >= 1.0e20) {
3727	columnUpperWork_[i] = COIN_DBL_MAX;
3728	} else {
3729	columnUpperWork_[i] = upperValue * multiplier;
3730	if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3731	if (columnLowerWork_[i] >= 0.0) {
3732	columnUpperWork_[i] = columnLowerWork_[i];
3733	} else if (columnUpperWork_[i] <= 0.0) {
3734	columnLowerWork_[i] = columnUpperWork_[i];
3735	} else {
3736	columnUpperWork_[i] = 0.0;
3737	columnLowerWork_[i] = 0.0;
3738	}
3739	}
3740	}
3741	} else if (upperValue < 1.0e20) {
3742	columnLowerWork_[i] = -COIN_DBL_MAX;
3743	columnUpperWork_[i] = upperValue * multiplier;
3744	} else {
3745	// free
3746	columnLowerWork_[i] = -COIN_DBL_MAX;
3747	columnUpperWork_[i] = COIN_DBL_MAX;
3748	}
3749	}
3750	}
3751	for (i = 0; i < numberRows_; i++) {
3752	double multiplier = rhsScale_ * rowScale[i];
3753	double lowerValue = rowLower_[i];
3754	double upperValue = rowUpper_[i];
3755	if (lowerValue > -1.0e20) {
3756	rowLowerWork_[i] = lowerValue * multiplier;
3757	if (upperValue >= 1.0e20) {
3758	rowUpperWork_[i] = COIN_DBL_MAX;
3759	} else {
3760	rowUpperWork_[i] = upperValue * multiplier;
3761	if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3762	if (rowLowerWork_[i] >= 0.0) {
3763	rowUpperWork_[i] = rowLowerWork_[i];
3764	} else if (rowUpperWork_[i] <= 0.0) {
3765	rowLowerWork_[i] = rowUpperWork_[i];
3766	} else {
3767	rowUpperWork_[i] = 0.0;
3768	rowLowerWork_[i] = 0.0;
3769	}
3770	}
3771	}
3772	} else if (upperValue < 1.0e20) {
3773	rowLowerWork_[i] = -COIN_DBL_MAX;
3774	rowUpperWork_[i] = upperValue * multiplier;
3775	} else {
3776	// free
3777	rowLowerWork_[i] = -COIN_DBL_MAX;
3778	rowUpperWork_[i] = COIN_DBL_MAX;
3779	}
3780	}
3781	} else if (rhsScale_ != 1.0) {
3782	for (i = 0; i < numberColumns_; i++) {
3783	double lowerValue = columnLower_[i];
3784	double upperValue = columnUpper_[i];
3785	if (lowerValue > -1.0e20) {
3786	columnLowerWork_[i] = lowerValue * rhsScale_;
3787	if (upperValue >= 1.0e20) {
3788	columnUpperWork_[i] = COIN_DBL_MAX;
3789	} else {
3790	columnUpperWork_[i] = upperValue * rhsScale_;
3791	if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3792	if (columnLowerWork_[i] >= 0.0) {
3793	columnUpperWork_[i] = columnLowerWork_[i];
3794	} else if (columnUpperWork_[i] <= 0.0) {
3795	columnLowerWork_[i] = columnUpperWork_[i];
3796	} else {
3797	columnUpperWork_[i] = 0.0;
3798	columnLowerWork_[i] = 0.0;
3799	}
3800	}
3801	}
3802	} else if (upperValue < 1.0e20) {
3803	columnLowerWork_[i] = -COIN_DBL_MAX;
3804	columnUpperWork_[i] = upperValue * rhsScale_;
3805	} else {
3806	// free
3807	columnLowerWork_[i] = -COIN_DBL_MAX;
3808	columnUpperWork_[i] = COIN_DBL_MAX;
3809	}
3810	}
3811	for (i = 0; i < numberRows_; i++) {
3812	double lowerValue = rowLower_[i];
3813	double upperValue = rowUpper_[i];
3814	if (lowerValue > -1.0e20) {
3815	rowLowerWork_[i] = lowerValue * rhsScale_;
3816	if (upperValue >= 1.0e20) {
3817	rowUpperWork_[i] = COIN_DBL_MAX;
3818	} else {
3819	rowUpperWork_[i] = upperValue * rhsScale_;
3820	if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3821	if (rowLowerWork_[i] >= 0.0) {
3822	rowUpperWork_[i] = rowLowerWork_[i];
3823	} else if (rowUpperWork_[i] <= 0.0) {
3824	rowLowerWork_[i] = rowUpperWork_[i];
3825	} else {
3826	rowUpperWork_[i] = 0.0;
3827	rowLowerWork_[i] = 0.0;
3828	}
3829	}
3830	}
3831	} else if (upperValue < 1.0e20) {
3832	rowLowerWork_[i] = -COIN_DBL_MAX;
3833	rowUpperWork_[i] = upperValue * rhsScale_;
3834	} else {
3835	// free
3836	rowLowerWork_[i] = -COIN_DBL_MAX;
3837	rowUpperWork_[i] = COIN_DBL_MAX;
3838	}
3839	}
3840	} else {
3841	for (i = 0; i < numberColumns_; i++) {
3842	double lowerValue = columnLower_[i];
3843	double upperValue = columnUpper_[i];
3844	if (lowerValue > -1.0e20) {
3845	columnLowerWork_[i] = lowerValue;
3846	if (upperValue >= 1.0e20) {
3847	columnUpperWork_[i] = COIN_DBL_MAX;
3848	} else {
3849	columnUpperWork_[i] = upperValue;
3850	if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3851	if (columnLowerWork_[i] >= 0.0) {
3852	columnUpperWork_[i] = columnLowerWork_[i];
3853	} else if (columnUpperWork_[i] <= 0.0) {
3854	columnLowerWork_[i] = columnUpperWork_[i];
3855	} else {
3856	columnUpperWork_[i] = 0.0;
3857	columnLowerWork_[i] = 0.0;
3858	}
3859	}
3860	}
3861	} else if (upperValue < 1.0e20) {
3862	columnLowerWork_[i] = -COIN_DBL_MAX;
3863	columnUpperWork_[i] = upperValue;
3864	} else {
3865	// free
3866	columnLowerWork_[i] = -COIN_DBL_MAX;
3867	columnUpperWork_[i] = COIN_DBL_MAX;
3868	}
3869	}
3870	for (i = 0; i < numberRows_; i++) {
3871	double lowerValue = rowLower_[i];
3872	double upperValue = rowUpper_[i];
3873	if (lowerValue > -1.0e20) {
3874	rowLowerWork_[i] = lowerValue;
3875	if (upperValue >= 1.0e20) {
3876	rowUpperWork_[i] = COIN_DBL_MAX;
3877	} else {
3878	rowUpperWork_[i] = upperValue;
3879	if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
3880	if (rowLowerWork_[i] >= 0.0) {
3881	rowUpperWork_[i] = rowLowerWork_[i];
3882	} else if (rowUpperWork_[i] <= 0.0) {
3883	rowLowerWork_[i] = rowUpperWork_[i];
3884	} else {
3885	rowUpperWork_[i] = 0.0;
3886	rowLowerWork_[i] = 0.0;
3887	}
3888	}
3889	}
3890	} else if (upperValue < 1.0e20) {
3891	rowLowerWork_[i] = -COIN_DBL_MAX;
3892	rowUpperWork_[i] = upperValue;
3893	} else {
3894	// free
3895	rowLowerWork_[i] = -COIN_DBL_MAX;
3896	rowUpperWork_[i] = COIN_DBL_MAX;
3897	}
3898	}
3899	}
3900	}
3901	if (what == 63) {
3902	// move information to work arrays
3903	double direction = optimizationDirection_;
3904	// direction is actually scale out not scale in
3905	if (direction)
3906	direction = 1.0 / direction;
3907	if (direction != 1.0) {
3908	// reverse all dual signs
3909	for (i = 0; i < numberColumns_; i++)
3910	reducedCost_[i] *= direction;
3911	for (i = 0; i < numberRows_; i++)
3912	dual_[i] *= direction;
3913	}
3914	for (i = 0; i < numberRows_ + numberColumns_; i++) {
3915	setFakeBound(i, noFake);
3916	}
3917	if (rowScale_) {
3918	const double * obj = objective();
3919	double direction = optimizationDirection_ * objectiveScale_;
3920	// clean up any mismatches on infinity
3921	// and fix any variables with tiny gaps
3922	double primalTolerance = dblParam_[ClpPrimalTolerance];
3923	// on entry
3924	const double * inverseScale = inverseColumnScale_;
3925	for (i = 0; i < numberColumns_; i++) {
3926	CoinAssert(fabs(obj[i]) < 1.0e25);
3927	double scaleFactor = columnScale_[i];
3928	double multiplier = rhsScale_ * inverseScale[i];
3929	scaleFactor *= direction;
3930	objectiveWork_[i] = obj[i] * scaleFactor;
3931	reducedCostWork_[i] = reducedCost_[i] * scaleFactor;
3932	double lowerValue = columnLower_[i];
3933	double upperValue = columnUpper_[i];
3934	if (lowerValue > -1.0e20) {
3935	columnLowerWork_[i] = lowerValue * multiplier;
3936	if (upperValue >= 1.0e20) {
3937	columnUpperWork_[i] = COIN_DBL_MAX;
3938	} else {
3939	columnUpperWork_[i] = upperValue * multiplier;
3940	if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
3941	if (columnLowerWork_[i] >= 0.0) {
3942	columnUpperWork_[i] = columnLowerWork_[i];
3943	} else if (columnUpperWork_[i] <= 0.0) {
3944	columnLowerWork_[i] = columnUpperWork_[i];
3945	} else {
3946	columnUpperWork_[i] = 0.0;
3947	columnLowerWork_[i] = 0.0;
3948	}
3949	}
3950	}
3951	} else if (upperValue < 1.0e20) {
3952	columnLowerWork_[i] = -COIN_DBL_MAX;
3953	columnUpperWork_[i] = upperValue * multiplier;
3954	} else {
3955	// free
3956	columnLowerWork_[i] = -COIN_DBL_MAX;
3957	columnUpperWork_[i] = COIN_DBL_MAX;
3958	}
3959	double value = columnActivity_[i] * multiplier;
3960	if (fabs(value) > 1.0e20) {
3961	//printf("bad value of %g for column %d\n",value,i);
3962	setColumnStatus(i, superBasic);
3963	if (columnUpperWork_[i] < 0.0) {
3964	value = columnUpperWork_[i];
3965	} else if (columnLowerWork_[i] > 0.0) {
3966	value = columnLowerWork_[i];
3967	} else {
3968	value = 0.0;
3969	}
3970	}
3971	columnActivityWork_[i] = value;
3972	}
3973	inverseScale = inverseRowScale_;
3974	for (i = 0; i < numberRows_; i++) {
3975	dual_[i] *= inverseScale[i];
3976	dual_[i] *= objectiveScale_;
3977	rowReducedCost_[i] = dual_[i];
3978	double multiplier = rhsScale_ * rowScale_[i];
3979	double value = rowActivity_[i] * multiplier;
3980	if (fabs(value) > 1.0e20) {
3981	//printf("bad value of %g for row %d\n",value,i);
3982	setRowStatus(i, superBasic);
3983	if (rowUpperWork_[i] < 0.0) {
3984	value = rowUpperWork_[i];
3985	} else if (rowLowerWork_[i] > 0.0) {
3986	value = rowLowerWork_[i];
3987	} else {
3988	value = 0.0;
3989	}
3990	}
3991	rowActivityWork_[i] = value;
3992	}
3993	} else if (objectiveScale_ != 1.0 \|\| rhsScale_ != 1.0) {
3994	// on entry
3995	for (i = 0; i < numberColumns_; i++) {
3996	double value = columnActivity_[i];
3997	value *= rhsScale_;
3998	if (fabs(value) > 1.0e20) {
3999	//printf("bad value of %g for column %d\n",value,i);
4000	setColumnStatus(i, superBasic);
4001	if (columnUpperWork_[i] < 0.0) {
4002	value = columnUpperWork_[i];
4003	} else if (columnLowerWork_[i] > 0.0) {
4004	value = columnLowerWork_[i];
4005	} else {
4006	value = 0.0;
4007	}
4008	}
4009	columnActivityWork_[i] = value;
4010	reducedCostWork_[i] = reducedCost_[i] * objectiveScale_;
4011	}
4012	for (i = 0; i < numberRows_; i++) {
4013	double value = rowActivity_[i];
4014	value *= rhsScale_;
4015	if (fabs(value) > 1.0e20) {
4016	//printf("bad value of %g for row %d\n",value,i);
4017	setRowStatus(i, superBasic);
4018	if (rowUpperWork_[i] < 0.0) {
4019	value = rowUpperWork_[i];
4020	} else if (rowLowerWork_[i] > 0.0) {
4021	value = rowLowerWork_[i];
4022	} else {
4023	value = 0.0;
4024	}
4025	}
4026	rowActivityWork_[i] = value;
4027	dual_[i] *= objectiveScale_;
4028	rowReducedCost_[i] = dual_[i];
4029	}
4030	} else {
4031	// on entry
4032	for (i = 0; i < numberColumns_; i++) {
4033	double value = columnActivity_[i];
4034	if (fabs(value) > 1.0e20) {
4035	//printf("bad value of %g for column %d\n",value,i);
4036	setColumnStatus(i, superBasic);
4037	if (columnUpperWork_[i] < 0.0) {
4038	value = columnUpperWork_[i];
4039	} else if (columnLowerWork_[i] > 0.0) {
4040	value = columnLowerWork_[i];
4041	} else {
4042	value = 0.0;
4043	}
4044	}
4045	columnActivityWork_[i] = value;
4046	reducedCostWork_[i] = reducedCost_[i];
4047	}
4048	for (i = 0; i < numberRows_; i++) {
4049	double value = rowActivity_[i];
4050	if (fabs(value) > 1.0e20) {
4051	//printf("bad value of %g for row %d\n",value,i);
4052	setRowStatus(i, superBasic);
4053	if (rowUpperWork_[i] < 0.0) {
4054	value = rowUpperWork_[i];
4055	} else if (rowLowerWork_[i] > 0.0) {
4056	value = rowLowerWork_[i];
4057	} else {
4058	value = 0.0;
4059	}
4060	}
4061	rowActivityWork_[i] = value;
4062	rowReducedCost_[i] = dual_[i];
4063	}
4064	}
4065	}
4066
4067	if (what == 63 && doSanityCheck) {
4068	// check rim of problem okay
4069	if (!sanityCheck())
4070	goodMatrix = false;
4071	}
4072	// we need to treat matrix as if each element by rowScaleIn and columnScaleout??
4073	// maybe we need to move scales to SimplexModel for factorization?
4074	if ((what == 63 && !pivotVariable_) \|\| (newArrays && !keepPivots)) {
4075	delete [] pivotVariable_;
4076	pivotVariable_ = new int[numberRows2];
4077	for (int i = 0; i < numberRows2; i++)
4078	pivotVariable_[i] = -1;
4079	} else if (what == 63 && !keepPivots) {
4080	// just reset
4081	for (int i = 0; i < numberRows2; i++)
4082	pivotVariable_[i] = -1;
4083	} else if (what == 63) {
4084	// check pivots
4085	for (int i = 0; i < numberRows2; i++) {
4086	int iSequence = pivotVariable_[i];
4087	if (iSequence<0\|\|(iSequence < numberRows_ + numberColumns_ &&
4088	getStatus(iSequence) != basic)) {
4089	keepPivots = false;
4090	break;
4091	}
4092	}
4093	if (!keepPivots) {
4094	// reset
4095	for (int i = 0; i < numberRows2; i++)
4096	pivotVariable_[i] = -1;
4097	} else {
4098	// clean
4099	for (int i = 0; i < numberColumns_ + numberRows_; i++) {
4100	Status status = getStatus(i);
4101	if (status != basic) {
4102	if (upper_[i] == lower_[i]) {
4103	setStatus(i, isFixed);
4104	solution_[i] = lower_[i];
4105	} else if (status == atLowerBound) {
4106	if (lower_[i] > -1.0e20) {
4107	solution_[i] = lower_[i];
4108	} else {
4109	//printf("seq %d at lower of %g\n",i,lower_[i]);
4110	if (upper_[i] < 1.0e20) {
4111	solution_[i] = upper_[i];
4112	setStatus(i, atUpperBound);
4113	} else {
4114	setStatus(i, isFree);
4115	}
4116	}
4117	} else if (status == atUpperBound) {
4118	if (upper_[i] < 1.0e20) {
4119	solution_[i] = upper_[i];
4120	} else {
4121	//printf("seq %d at upper of %g\n",i,upper_[i]);
4122	if (lower_[i] > -1.0e20) {
4123	solution_[i] = lower_[i];
4124	setStatus(i, atLowerBound);
4125	} else {
4126	setStatus(i, isFree);
4127	}
4128	}
4129	} else if (status == isFixed && upper_[i] > lower_[i]) {
4130	// was fixed - not now
4131	if (solution_[i] <= lower_[i]+primalTolerance_) {
4132	setStatus(i, atLowerBound);
4133	} else if (solution_[i] >= upper_[i]-primalTolerance_) {
4134	setStatus(i, atUpperBound);
4135	} else {
4136	setStatus(i, superBasic);
4137	}
4138	}
4139	}
4140	}
4141	}
4142	}
4143
4144	if (what == 63) {
4145	if (newArrays) {
4146	// get some arrays
4147	int iRow, iColumn;
4148	// these are "indexed" arrays so we always know where nonzeros are
4149	/**********************************************************
4150	rowArray_[3] is long enough for rows+columns
4151	rowArray_[1] is long enough for max(rows,columns)
4152	*********************************************************/
4153	for (iRow = 0; iRow < 4; iRow++) {
4154	int length = numberRows2 + factorization_->maximumPivots();
4155	if (iRow == 3 \|\| objective_->type() > 1)
4156	length += numberColumns_;
4157	else if (iRow == 1)
4158	length = CoinMax(length, numberColumns_);
4159	if ((specialOptions_ & 65536) == 0 \|\| !rowArray_[iRow]) {
4160	delete rowArray_[iRow];
4161	rowArray_[iRow] = new CoinIndexedVector();
4162	}
4163	rowArray_[iRow]->reserve(length);
4164	}
4165
4166	for (iColumn = 0; iColumn < 2; iColumn++) {
4167	if ((specialOptions_ & 65536) == 0 \|\| !columnArray_[iColumn]) {
4168	delete columnArray_[iColumn];
4169	columnArray_[iColumn] = new CoinIndexedVector();
4170	}
4171	columnArray_[iColumn]->reserve(numberColumns_+numberRows2);
4172	}
4173	} else {
4174	int iRow, iColumn;
4175	for (iRow = 0; iRow < 4; iRow++) {
4176	int length = numberRows2 + factorization_->maximumPivots();
4177	if (iRow == 3 \|\| objective_->type() > 1)
4178	length += numberColumns_;
4179	if(rowArray_[iRow]->capacity() >= length) {
4180	rowArray_[iRow]->clear();
4181	} else {
4182	// model size or maxinv changed
4183	rowArray_[iRow]->reserve(length);
4184	}
4185	#ifndef NDEBUG
4186	rowArray_[iRow]->checkClear();
4187	#endif
4188	}
4189
4190	for (iColumn = 0; iColumn < 2; iColumn++) {
4191	int length = numberColumns_;
4192	if (iColumn)
4193	length = CoinMax(numberRows2, numberColumns_);
4194	if(columnArray_[iColumn]->capacity() >= length) {
4195	columnArray_[iColumn]->clear();
4196	} else {
4197	// model size or maxinv changed
4198	columnArray_[iColumn]->reserve(length);
4199	}
4200	#ifndef NDEBUG
4201	columnArray_[iColumn]->checkClear();
4202	#endif
4203	}
4204	}
4205	}
4206	if (problemStatus_ == 10) {
4207	problemStatus_ = -1;
4208	handler_->setLogLevel(saveLevel); // switch back messages
4209	}
4210	if ((what & 5) != 0)
4211	matrix_->generalExpanded(this, 9, what); // update costs and bounds if necessary
4212	if (goodMatrix && (specialOptions_ & 65536) != 0) {
4213	int save = maximumColumns_ + maximumRows_;
4214	CoinMemcpyN(cost_, numberTotal, cost_ + save);
4215	CoinMemcpyN(lower_, numberTotal, lower_ + save);
4216	CoinMemcpyN(upper_, numberTotal, upper_ + save);
4217	CoinMemcpyN(dj_, numberTotal, dj_ + save);
4218	CoinMemcpyN(solution_, numberTotal, solution_ + save);
4219	if (rowScale_ && !savedRowScale_) {
4220	double * temp;
4221	temp = new double [4*maximumRows_];
4222	CoinFillN(temp, 4 * maximumRows_, 1.0);
4223	CoinMemcpyN(rowScale_, numberRows2, temp);
4224	CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + maximumRows_);
4225	CoinMemcpyN(rowScale_, numberRows2, temp + 2 * maximumRows_);
4226	CoinMemcpyN(rowScale_ + numberRows2, numberRows2, temp + 3 * maximumRows_);
4227	delete [] rowScale_;
4228	savedRowScale_ = temp;
4229	rowScale_ = savedRowScale_;
4230	inverseRowScale_ = savedRowScale_ + maximumInternalRows_;
4231	temp = new double [4*maximumColumns_];
4232	CoinFillN(temp, 4 * maximumColumns_, 1.0);
4233	CoinMemcpyN(columnScale_, numberColumns_, temp);
4234	CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + maximumColumns_);
4235	CoinMemcpyN(columnScale_, numberColumns_, temp + 2 * maximumColumns_);
4236	CoinMemcpyN(columnScale_ + numberColumns_, numberColumns_, temp + 3 * maximumColumns_);
4237	delete [] columnScale_;
4238	savedColumnScale_ = temp;
4239	columnScale_ = savedColumnScale_;
4240	inverseColumnScale_ = savedColumnScale_ + maximumInternalColumns_;
4241	}
4242	}
4243	#ifdef CLP_USER_DRIVEN
4244	eventHandler_->event(ClpEventHandler::endOfCreateRim);
4245	#endif
4246	return goodMatrix;
4247	}
4248	// Does rows and columns
4249	void
4250	ClpSimplex::createRim1(bool initial)
4251	{
4252	int i;
4253	int numberRows2 = numberRows_ + numberExtraRows_;
4254	int numberTotal = numberRows2 + numberColumns_;
4255	if ((specialOptions_ & 65536) != 0 && true) {
4256	assert (!initial);
4257	int save = maximumColumns_ + maximumRows_;
4258	CoinMemcpyN(lower_ + save, numberTotal, lower_);
4259	CoinMemcpyN(upper_ + save, numberTotal, upper_);
4260	return;
4261	}
4262	const double * rowScale = rowScale_;
4263	// clean up any mismatches on infinity
4264	// and fix any variables with tiny gaps
4265	double primalTolerance = dblParam_[ClpPrimalTolerance];
4266	if(rowScale) {
4267	// If scaled then do all columns later in one loop
4268	if (!initial) {
4269	const double * inverseScale = inverseColumnScale_;
4270	for (i = 0; i < numberColumns_; i++) {
4271	double multiplier = rhsScale_ * inverseScale[i];
4272	double lowerValue = columnLower_[i];
4273	double upperValue = columnUpper_[i];
4274	if (lowerValue > -1.0e20) {
4275	columnLowerWork_[i] = lowerValue * multiplier;
4276	if (upperValue >= 1.0e20) {
4277	columnUpperWork_[i] = COIN_DBL_MAX;
4278	} else {
4279	columnUpperWork_[i] = upperValue * multiplier;
4280	if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4281	if (columnLowerWork_[i] >= 0.0) {
4282	columnUpperWork_[i] = columnLowerWork_[i];
4283	} else if (columnUpperWork_[i] <= 0.0) {
4284	columnLowerWork_[i] = columnUpperWork_[i];
4285	} else {
4286	columnUpperWork_[i] = 0.0;
4287	columnLowerWork_[i] = 0.0;
4288	}
4289	}
4290	}
4291	} else if (upperValue < 1.0e20) {
4292	columnLowerWork_[i] = -COIN_DBL_MAX;
4293	columnUpperWork_[i] = upperValue * multiplier;
4294	} else {
4295	// free
4296	columnLowerWork_[i] = -COIN_DBL_MAX;
4297	columnUpperWork_[i] = COIN_DBL_MAX;
4298	}
4299	}
4300	}
4301	for (i = 0; i < numberRows_; i++) {
4302	double multiplier = rhsScale_ * rowScale[i];
4303	double lowerValue = rowLower_[i];
4304	double upperValue = rowUpper_[i];
4305	if (lowerValue > -1.0e20) {
4306	rowLowerWork_[i] = lowerValue * multiplier;
4307	if (upperValue >= 1.0e20) {
4308	rowUpperWork_[i] = COIN_DBL_MAX;
4309	} else {
4310	rowUpperWork_[i] = upperValue * multiplier;
4311	if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4312	if (rowLowerWork_[i] >= 0.0) {
4313	rowUpperWork_[i] = rowLowerWork_[i];
4314	} else if (rowUpperWork_[i] <= 0.0) {
4315	rowLowerWork_[i] = rowUpperWork_[i];
4316	} else {
4317	rowUpperWork_[i] = 0.0;
4318	rowLowerWork_[i] = 0.0;
4319	}
4320	}
4321	}
4322	} else if (upperValue < 1.0e20) {
4323	rowLowerWork_[i] = -COIN_DBL_MAX;
4324	rowUpperWork_[i] = upperValue * multiplier;
4325	} else {
4326	// free
4327	rowLowerWork_[i] = -COIN_DBL_MAX;
4328	rowUpperWork_[i] = COIN_DBL_MAX;
4329	}
4330	}
4331	} else if (rhsScale_ != 1.0) {
4332	for (i = 0; i < numberColumns_; i++) {
4333	double lowerValue = columnLower_[i];
4334	double upperValue = columnUpper_[i];
4335	if (lowerValue > -1.0e20) {
4336	columnLowerWork_[i] = lowerValue * rhsScale_;
4337	if (upperValue >= 1.0e20) {
4338	columnUpperWork_[i] = COIN_DBL_MAX;
4339	} else {
4340	columnUpperWork_[i] = upperValue * rhsScale_;
4341	if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4342	if (columnLowerWork_[i] >= 0.0) {
4343	columnUpperWork_[i] = columnLowerWork_[i];
4344	} else if (columnUpperWork_[i] <= 0.0) {
4345	columnLowerWork_[i] = columnUpperWork_[i];
4346	} else {
4347	columnUpperWork_[i] = 0.0;
4348	columnLowerWork_[i] = 0.0;
4349	}
4350	}
4351	}
4352	} else if (upperValue < 1.0e20) {
4353	columnLowerWork_[i] = -COIN_DBL_MAX;
4354	columnUpperWork_[i] = upperValue * rhsScale_;
4355	} else {
4356	// free
4357	columnLowerWork_[i] = -COIN_DBL_MAX;
4358	columnUpperWork_[i] = COIN_DBL_MAX;
4359	}
4360	}
4361	for (i = 0; i < numberRows_; i++) {
4362	double lowerValue = rowLower_[i];
4363	double upperValue = rowUpper_[i];
4364	if (lowerValue > -1.0e20) {
4365	rowLowerWork_[i] = lowerValue * rhsScale_;
4366	if (upperValue >= 1.0e20) {
4367	rowUpperWork_[i] = COIN_DBL_MAX;
4368	} else {
4369	rowUpperWork_[i] = upperValue * rhsScale_;
4370	if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4371	if (rowLowerWork_[i] >= 0.0) {
4372	rowUpperWork_[i] = rowLowerWork_[i];
4373	} else if (rowUpperWork_[i] <= 0.0) {
4374	rowLowerWork_[i] = rowUpperWork_[i];
4375	} else {
4376	rowUpperWork_[i] = 0.0;
4377	rowLowerWork_[i] = 0.0;
4378	}
4379	}
4380	}
4381	} else if (upperValue < 1.0e20) {
4382	rowLowerWork_[i] = -COIN_DBL_MAX;
4383	rowUpperWork_[i] = upperValue * rhsScale_;
4384	} else {
4385	// free
4386	rowLowerWork_[i] = -COIN_DBL_MAX;
4387	rowUpperWork_[i] = COIN_DBL_MAX;
4388	}
4389	}
4390	} else {
4391	for (i = 0; i < numberColumns_; i++) {
4392	double lowerValue = columnLower_[i];
4393	double upperValue = columnUpper_[i];
4394	if (lowerValue > -1.0e20) {
4395	columnLowerWork_[i] = lowerValue;
4396	if (upperValue >= 1.0e20) {
4397	columnUpperWork_[i] = COIN_DBL_MAX;
4398	} else {
4399	columnUpperWork_[i] = upperValue;
4400	if (fabs(columnUpperWork_[i] - columnLowerWork_[i]) <= primalTolerance) {
4401	if (columnLowerWork_[i] >= 0.0) {
4402	columnUpperWork_[i] = columnLowerWork_[i];
4403	} else if (columnUpperWork_[i] <= 0.0) {
4404	columnLowerWork_[i] = columnUpperWork_[i];
4405	} else {
4406	columnUpperWork_[i] = 0.0;
4407	columnLowerWork_[i] = 0.0;
4408	}
4409	}
4410	}
4411	} else if (upperValue < 1.0e20) {
4412	columnLowerWork_[i] = -COIN_DBL_MAX;
4413	columnUpperWork_[i] = upperValue;
4414	} else {
4415	// free
4416	columnLowerWork_[i] = -COIN_DBL_MAX;
4417	columnUpperWork_[i] = COIN_DBL_MAX;
4418	}
4419	}
4420	for (i = 0; i < numberRows_; i++) {
4421	double lowerValue = rowLower_[i];
4422	double upperValue = rowUpper_[i];
4423	if (lowerValue > -1.0e20) {
4424	rowLowerWork_[i] = lowerValue;
4425	if (upperValue >= 1.0e20) {
4426	rowUpperWork_[i] = COIN_DBL_MAX;
4427	} else {
4428	rowUpperWork_[i] = upperValue;
4429	if (fabs(rowUpperWork_[i] - rowLowerWork_[i]) <= primalTolerance) {
4430	if (rowLowerWork_[i] >= 0.0) {
4431	rowUpperWork_[i] = rowLowerWork_[i];
4432	} else if (rowUpperWork_[i] <= 0.0) {
4433	rowLowerWork_[i] = rowUpperWork_[i];
4434	} else {
4435	rowUpperWork_[i] = 0.0;
4436	rowLowerWork_[i] = 0.0;
4437	}
4438	}
4439	}
4440	} else if (upperValue < 1.0e20) {
4441	rowLowerWork_[i] = -COIN_DBL_MAX;
4442	rowUpperWork_[i] = upperValue;
4443	} else {
4444	// free
4445	rowLowerWork_[i] = -COIN_DBL_MAX;
4446	rowUpperWork_[i] = COIN_DBL_MAX;
4447	}
4448	}
4449	}
4450	#ifndef NDEBUG
4451	if ((specialOptions_ & 65536) != 0 && false) {
4452	assert (!initial);
4453	int save = maximumColumns_ + maximumRows_;
4454	for (int i = 0; i < numberTotal; i++) {
4455	assert (fabs(lower_[i] - lower_[i+save]) < 1.0e-5);
4456	assert (fabs(upper_[i] - upper_[i+save]) < 1.0e-5);
4457	}
4458	}
4459	#endif
4460	}
4461	// Does objective
4462	void
4463	ClpSimplex::createRim4(bool initial)
4464	{
4465	int i;
4466	int numberRows2 = numberRows_ + numberExtraRows_;
4467	int numberTotal = numberRows2 + numberColumns_;
4468	if ((specialOptions_ & 65536) != 0 && true) {
4469	assert (!initial);
4470	int save = maximumColumns_ + maximumRows_;
4471	CoinMemcpyN(cost_ + save, numberTotal, cost_);
4472	return;
4473	}
4474	double direction = optimizationDirection_ * objectiveScale_;
4475	const double * obj = objective();
4476	const double * rowScale = rowScale_;
4477	const double * columnScale = columnScale_;
4478	// and also scale by scale factors
4479	if (rowScale) {
4480	if (rowObjective_) {
4481	for (i = 0; i < numberRows_; i++)
4482	rowObjectiveWork_[i] = rowObjective_[i] * direction / rowScale[i];
4483	} else {
4484	memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4485	}
4486	// If scaled then do all columns later in one loop
4487	if (!initial) {
4488	for (i = 0; i < numberColumns_; i++) {
4489	CoinAssert(fabs(obj[i]) < 1.0e25);
4490	objectiveWork_[i] = obj[i] * direction * columnScale[i];
4491	}
4492	}
4493	} else {
4494	if (rowObjective_) {
4495	for (i = 0; i < numberRows_; i++)
4496	rowObjectiveWork_[i] = rowObjective_[i] * direction;
4497	} else {
4498	memset(rowObjectiveWork_, 0, numberRows_ * sizeof(double));
4499	}
4500	for (i = 0; i < numberColumns_; i++) {
4501	CoinAssert(fabs(obj[i]) < 1.0e25);
4502	objectiveWork_[i] = obj[i] * direction;
4503	}
4504	}
4505	}
4506	// Does rows and columns and objective
4507	void
4508	ClpSimplex::createRim5(bool initial)
4509	{
4510	createRim4(initial);
4511	createRim1(initial);
4512	}
4513	void
4514	ClpSimplex::deleteRim(int getRidOfFactorizationData)
4515	{
4516	// Just possible empty problem
4517	int numberRows = numberRows_;
4518	int numberColumns = numberColumns_;
4519	if (!numberRows \|\| !numberColumns) {
4520	numberRows = 0;
4521	if (objective_->type() < 2)
4522	numberColumns = 0;
4523	}
4524	int i;
4525	if (problemStatus_ != 1 && problemStatus_ != 2) {
4526	delete [] ray_;
4527	ray_ = NULL;
4528	}
4529	// set upperOut_ to furthest away from bound so can use in dual for dualBound_
4530	upperOut_ = 1.0;
4531	#if 0
4532	{
4533	int nBad = 0;
4534	for (i = 0; i < numberColumns; i++) {
4535	if (lower_[i] == upper_[i] && getColumnStatus(i) == basic)
4536	nBad++;
4537	}
4538	if (nBad)
4539	printf("yy %d basic fixed\n", nBad);
4540	}
4541	#endif
4542	if ((moreSpecialOptions_&4194304)!=0) {
4543	// preset tolerances were changed
4544	moreSpecialOptions_ &= ~4194304;
4545	primalTolerance_=1.0e-7;
4546	dblParam_[ClpPrimalTolerance]=primalTolerance_;
4547	dualTolerance_=1.0e-7;
4548	dblParam_[ClpDualTolerance]=dualTolerance_;
4549	}
4550	// ray may be null if in branch and bound
4551	if (rowScale_) {
4552	// Collect infeasibilities
4553	int numberPrimalScaled = 0;
4554	int numberPrimalUnscaled = 0;
4555	int numberDualScaled = 0;
4556	int numberDualUnscaled = 0;
4557	double scaleC = 1.0 / objectiveScale_;
4558	double scaleR = 1.0 / rhsScale_;
4559	const double * inverseScale = inverseColumnScale_;
4560	for (i = 0; i < numberColumns; i++) {
4561	double scaleFactor = columnScale_[i];
4562	double valueScaled = columnActivityWork_[i];
4563	double lowerScaled = columnLowerWork_[i];
4564	double upperScaled = columnUpperWork_[i];
4565	if (lowerScaled > -1.0e20 \|\| upperScaled < 1.0e20) {
4566	if (valueScaled < lowerScaled - primalTolerance_ \|\|
4567	valueScaled > upperScaled + primalTolerance_)
4568	numberPrimalScaled++;
4569	else
4570	upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4571	}
4572	columnActivity_[i] = valueScaled * scaleFactor * scaleR;
4573	double value = columnActivity_[i];
4574	if (value < columnLower_[i] - primalTolerance_)
4575	numberPrimalUnscaled++;
4576	else if (value > columnUpper_[i] + primalTolerance_)
4577	numberPrimalUnscaled++;
4578	double valueScaledDual = reducedCostWork_[i];
4579	if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4580	numberDualScaled++;
4581	if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4582	numberDualScaled++;
4583	reducedCost_[i] = (valueScaledDual * scaleC) * inverseScale[i];
4584	double valueDual = reducedCost_[i];
4585	if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4586	numberDualUnscaled++;
4587	if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4588	numberDualUnscaled++;
4589	}
4590	inverseScale = inverseRowScale_;
4591	for (i = 0; i < numberRows; i++) {
4592	double scaleFactor = rowScale_[i];
4593	double valueScaled = rowActivityWork_[i];
4594	double lowerScaled = rowLowerWork_[i];
4595	double upperScaled = rowUpperWork_[i];
4596	if (lowerScaled > -1.0e20 \|\| upperScaled < 1.0e20) {
4597	if (valueScaled < lowerScaled - primalTolerance_ \|\|
4598	valueScaled > upperScaled + primalTolerance_)
4599	numberPrimalScaled++;
4600	else
4601	upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4602	}
4603	rowActivity_[i] = (valueScaled * scaleR) * inverseScale[i];
4604	double value = rowActivity_[i];
4605	if (value < rowLower_[i] - primalTolerance_)
4606	numberPrimalUnscaled++;
4607	else if (value > rowUpper_[i] + primalTolerance_)
4608	numberPrimalUnscaled++;
4609	double valueScaledDual = dual_[i] + rowObjectiveWork_[i];;
4610	if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4611	numberDualScaled++;
4612	if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4613	numberDualScaled++;
4614	dual_[i] = scaleFactor scaleC;
4615	double valueDual = dual_[i];
4616	if (rowObjective_)
4617	valueDual += rowObjective_[i];
4618	if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4619	numberDualUnscaled++;
4620	if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4621	numberDualUnscaled++;
4622	}
4623	if (!problemStatus_ && !secondaryStatus_) {
4624	// See if we need to set secondary status
4625	if (numberPrimalUnscaled) {
4626	if (numberDualUnscaled)
4627	secondaryStatus_ = 4;
4628	else
4629	secondaryStatus_ = 2;
4630	} else {
4631	if (numberDualUnscaled)
4632	secondaryStatus_ = 3;
4633	}
4634	}
4635	if (problemStatus_ == 2 && ray_) {
4636	for (i = 0; i < numberColumns; i++) {
4637	ray_[i] *= columnScale_[i];
4638	}
4639	} else if (problemStatus_ == 1 && ray_) {
4640	for (i = 0; i < numberRows; i++) {
4641	ray_[i] *= rowScale_[i];
4642	}
4643	}
4644	} else if (rhsScale_ != 1.0 \|\| objectiveScale_ != 1.0) {
4645	// Collect infeasibilities
4646	int numberPrimalScaled = 0;
4647	int numberPrimalUnscaled = 0;
4648	int numberDualScaled = 0;
4649	int numberDualUnscaled = 0;
4650	double scaleC = 1.0 / objectiveScale_;
4651	double scaleR = 1.0 / rhsScale_;
4652	for (i = 0; i < numberColumns; i++) {
4653	double valueScaled = columnActivityWork_[i];
4654	double lowerScaled = columnLowerWork_[i];
4655	double upperScaled = columnUpperWork_[i];
4656	if (lowerScaled > -1.0e20 \|\| upperScaled < 1.0e20) {
4657	if (valueScaled < lowerScaled - primalTolerance_ \|\|
4658	valueScaled > upperScaled + primalTolerance_)
4659	numberPrimalScaled++;
4660	else
4661	upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4662	}
4663	columnActivity_[i] = valueScaled * scaleR;
4664	double value = columnActivity_[i];
4665	if (value < columnLower_[i] - primalTolerance_)
4666	numberPrimalUnscaled++;
4667	else if (value > columnUpper_[i] + primalTolerance_)
4668	numberPrimalUnscaled++;
4669	double valueScaledDual = reducedCostWork_[i];
4670	if (valueScaled > columnLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4671	numberDualScaled++;
4672	if (valueScaled < columnUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4673	numberDualScaled++;
4674	reducedCost_[i] = valueScaledDual * scaleC;
4675	double valueDual = reducedCost_[i];
4676	if (value > columnLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4677	numberDualUnscaled++;
4678	if (value < columnUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4679	numberDualUnscaled++;
4680	}
4681	for (i = 0; i < numberRows; i++) {
4682	double valueScaled = rowActivityWork_[i];
4683	double lowerScaled = rowLowerWork_[i];
4684	double upperScaled = rowUpperWork_[i];
4685	if (lowerScaled > -1.0e20 \|\| upperScaled < 1.0e20) {
4686	if (valueScaled < lowerScaled - primalTolerance_ \|\|
4687	valueScaled > upperScaled + primalTolerance_)
4688	numberPrimalScaled++;
4689	else
4690	upperOut_ = CoinMax(upperOut_, CoinMin(valueScaled - lowerScaled, upperScaled - valueScaled));
4691	}
4692	rowActivity_[i] = valueScaled * scaleR;
4693	double value = rowActivity_[i];
4694	if (value < rowLower_[i] - primalTolerance_)
4695	numberPrimalUnscaled++;
4696	else if (value > rowUpper_[i] + primalTolerance_)
4697	numberPrimalUnscaled++;
4698	double valueScaledDual = dual_[i] + rowObjectiveWork_[i];;
4699	if (valueScaled > rowLowerWork_[i] + primalTolerance_ && valueScaledDual > dualTolerance_)
4700	numberDualScaled++;
4701	if (valueScaled < rowUpperWork_[i] - primalTolerance_ && valueScaledDual < -dualTolerance_)
4702	numberDualScaled++;
4703	dual_[i] *= scaleC;
4704	double valueDual = dual_[i];
4705	if (rowObjective_)
4706	valueDual += rowObjective_[i];
4707	if (value > rowLower_[i] + primalTolerance_ && valueDual > dualTolerance_)
4708	numberDualUnscaled++;
4709	if (value < rowUpper_[i] - primalTolerance_ && valueDual < -dualTolerance_)
4710	numberDualUnscaled++;
4711	}
4712	if (!problemStatus_ && !secondaryStatus_) {
4713	// See if we need to set secondary status
4714	if (numberPrimalUnscaled) {
4715	if (numberDualUnscaled)
4716	secondaryStatus_ = 4;
4717	else
4718	secondaryStatus_ = 2;
4719	} else {
4720	if (numberDualUnscaled)
4721	secondaryStatus_ = 3;
4722	}
4723	}
4724	} else {
4725	if (columnActivityWork_) {
4726	for (i = 0; i < numberColumns; i++) {
4727	double value = columnActivityWork_[i];
4728	double lower = columnLowerWork_[i];
4729	double upper = columnUpperWork_[i];
4730	if (lower > -1.0e20 \|\| upper < 1.0e20) {
4731	if (value > lower && value < upper)
4732	upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value));
4733	}
4734	columnActivity_[i] = columnActivityWork_[i];
4735	reducedCost_[i] = reducedCostWork_[i];
4736	}
4737	for (i = 0; i < numberRows; i++) {
4738	double value = rowActivityWork_[i];
4739	double lower = rowLowerWork_[i];
4740	double upper = rowUpperWork_[i];
4741	if (lower > -1.0e20 \|\| upper < 1.0e20) {
4742	if (value > lower && value < upper)
4743	upperOut_ = CoinMax(upperOut_, CoinMin(value - lower, upper - value));
4744	}
4745	rowActivity_[i] = rowActivityWork_[i];
4746	}
4747	}
4748	}
4749	// switch off scalefactor if auto
4750	if (automaticScale_) {
4751	rhsScale_ = 1.0;
4752	objectiveScale_ = 1.0;
4753	}
4754	if (optimizationDirection_ != 1.0) {
4755	// and modify all dual signs
4756	for (i = 0; i < numberColumns; i++)
4757	reducedCost_[i] *= optimizationDirection_;
4758	for (i = 0; i < numberRows; i++)
4759	dual_[i] *= optimizationDirection_;
4760	}
4761	// scaling may have been turned off
4762	scalingFlag_ = abs(scalingFlag_);
4763	if(getRidOfFactorizationData > 0) {
4764	gutsOfDelete(getRidOfFactorizationData + 1);
4765	} else {
4766	// at least get rid of nonLinearCost_
4767	delete nonLinearCost_;
4768	nonLinearCost_ = NULL;
4769	}
4770	if (!rowObjective_ && problemStatus_ == 0 && objective_->type() == 1 &&
4771	numberRows && numberColumns) {
4772	// Redo objective value
4773	double objectiveValue = 0.0;
4774	const double * cost = objective();
4775	for (int i = 0; i < numberColumns; i++) {
4776	double value = columnActivity_[i];
4777	objectiveValue += value * cost[i];
4778	}
4779	//if (fabs(objectiveValue_ -objectiveValue*optimizationDirection())>1.0e-5)
4780	//printf("old obj %g new %g\n",objectiveValue_, objectiveValue*optimizationDirection());
4781	objectiveValue_ = objectiveValue * optimizationDirection();
4782	}
4783	// get rid of data
4784	matrix_->generalExpanded(this, 13, scalingFlag_);
4785	}
4786	void
4787	ClpSimplex::setDualBound(double value)
4788	{
4789	if (value > 0.0)
4790	dualBound_ = value;
4791	}
4792	void
4793	ClpSimplex::setInfeasibilityCost(double value)
4794	{
4795	if (value > 0.0)
4796	infeasibilityCost_ = value;
4797	}
4798	void ClpSimplex::setNumberRefinements( int value)
4799	{
4800	if (value >= 0 && value < 10)
4801	numberRefinements_ = value;
4802	}
4803	// Sets row pivot choice algorithm in dual
4804	void
4805	ClpSimplex::setDualRowPivotAlgorithm(ClpDualRowPivot & choice)
4806	{
4807	delete dualRowPivot_;
4808	dualRowPivot_ = choice.clone(true);
4809	dualRowPivot_->setModel(this);
4810	}
4811	// Sets column pivot choice algorithm in primal
4812	void
4813	ClpSimplex::setPrimalColumnPivotAlgorithm(ClpPrimalColumnPivot & choice)
4814	{
4815	delete primalColumnPivot_;
4816	primalColumnPivot_ = choice.clone(true);
4817	primalColumnPivot_->setModel(this);
4818	}
4819	void
4820	ClpSimplex::setFactorization( ClpFactorization & factorization)
4821	{
4822	if (factorization_)
4823	factorization_->setFactorization(factorization);
4824	else
4825	factorization_ = new ClpFactorization(factorization,
4826	numberRows_);
4827	}
4828
4829	// Swaps factorization
4830	ClpFactorization *
4831	ClpSimplex::swapFactorization( ClpFactorization * factorization)
4832	{
4833	ClpFactorization * swap = factorization_;
4834	factorization_ = factorization;
4835	return swap;
4836	}
4837	// Copies in factorization to existing one
4838	void
4839	ClpSimplex::copyFactorization( ClpFactorization & factorization)
4840	{
4841	*factorization_ = factorization;
4842	}
4843	/* Perturbation:
4844	-50 to +50 - perturb by this power of ten (-6 sounds good)
4845	100 - auto perturb if takes too long (1.0e-6 largest nonzero)
4846	101 - we are perturbed
4847	102 - don't try perturbing again
4848	default is 100
4849	*/
4850	void
4851	ClpSimplex::setPerturbation(int value)
4852	{
4853	if(value <= 100 && value >= -1000) {
4854	perturbation_ = value;
4855	}
4856	}
4857	// Sparsity on or off
4858	bool
4859	ClpSimplex::sparseFactorization() const
4860	{
4861	return factorization_->sparseThreshold() != 0;
4862	}
4863	void
4864	ClpSimplex::setSparseFactorization(bool value)
4865	{
4866	if (value) {
4867	if (!factorization_->sparseThreshold())
4868	factorization_->goSparse();
4869	} else {
4870	factorization_->sparseThreshold(0);
4871	}
4872	}
4873	void checkCorrect(ClpSimplex * /model/, int iRow,
4874	const double * element, const int * rowStart, const int * rowLength,
4875	const int * column,
4876	const double * columnLower_, const double * columnUpper_,
4877	int /infiniteUpperC/,
4878	int /infiniteLowerC/,
4879	double &maximumUpC,
4880	double &maximumDownC)
4881	{
4882	int infiniteUpper = 0;
4883	int infiniteLower = 0;
4884	double maximumUp = 0.0;
4885	double maximumDown = 0.0;
4886	CoinBigIndex rStart = rowStart[iRow];
4887	CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
4888	CoinBigIndex j;
4889	double large = 1.0e15;
4890	int iColumn;
4891	// Compute possible lower and upper ranges
4892
4893	for (j = rStart; j < rEnd; ++j) {
4894	double value = element[j];
4895	iColumn = column[j];
4896	if (value > 0.0) {
4897	if (columnUpper_[iColumn] >= large) {
4898	++infiniteUpper;
4899	} else {
4900	maximumUp += columnUpper_[iColumn] * value;
4901	}
4902	if (columnLower_[iColumn] <= -large) {
4903	++infiniteLower;
4904	} else {
4905	maximumDown += columnLower_[iColumn] * value;
4906	}
4907	} else if (value < 0.0) {
4908	if (columnUpper_[iColumn] >= large) {
4909	++infiniteLower;
4910	} else {
4911	maximumDown += columnUpper_[iColumn] * value;
4912	}
4913	if (columnLower_[iColumn] <= -large) {
4914	++infiniteUpper;
4915	} else {
4916	maximumUp += columnLower_[iColumn] * value;
4917	}
4918	}
4919	}
4920	//assert (infiniteLowerC==infiniteLower);
4921	//assert (infiniteUpperC==infiniteUpper);
4922	if (fabs(maximumUp - maximumUpC) > 1.0e-12 * CoinMax(fabs(maximumUp), fabs(maximumUpC)))
4923	COIN_DETAIL_PRINT(printf("row %d comp up %g, true up %g\n", iRow,
4924	maximumUpC, maximumUp));
4925	if (fabs(maximumDown - maximumDownC) > 1.0e-12 * CoinMax(fabs(maximumDown), fabs(maximumDownC)))
4926	COIN_DETAIL_PRINT(printf("row %d comp down %g, true down %g\n", iRow,
4927	maximumDownC, maximumDown));
4928	maximumUpC = maximumUp;
4929	maximumDownC = maximumDown;
4930	}
4931
4932	/* Tightens primal bounds to make dual faster. Unless
4933	fixed, bounds are slightly looser than they could be.
4934	This is to make dual go faster and is probably not needed
4935	with a presolve. Returns non-zero if problem infeasible
4936
4937	Fudge for branch and bound - put bounds on columns of factor *
4938	largest value (at continuous) - should improve stability
4939	in branch and bound on infeasible branches (0.0 is off)
4940	*/
4941	int
4942	ClpSimplex::tightenPrimalBounds(double factor, int doTight, bool tightIntegers)
4943	{
4944
4945	// Get a row copy in standard format
4946	CoinPackedMatrix copy;
4947	copy.setExtraGap(0.0);
4948	copy.setExtraMajor(0.0);
4949	copy.reverseOrderedCopyOf(*matrix());
4950	// Matrix may have been created so get rid of it
4951	matrix_->releasePackedMatrix();
4952	// get matrix data pointers
4953	const int * column = copy.getIndices();
4954	const CoinBigIndex * rowStart = copy.getVectorStarts();
4955	const int * rowLength = copy.getVectorLengths();
4956	const double * element = copy.getElements();
4957	int numberChanged = 1, iPass = 0;
4958	double large = largeValue(); // treat bounds > this as infinite
4959	#ifndef NDEBUG
4960	double large2 = 1.0e10 * large;
4961	#endif
4962	int numberInfeasible = 0;
4963	int totalTightened = 0;
4964
4965	double tolerance = primalTolerance();
4966
4967
4968	// Save column bounds
4969	double * saveLower = new double [numberColumns_];
4970	CoinMemcpyN(columnLower_, numberColumns_, saveLower);
4971	double * saveUpper = new double [numberColumns_];
4972	CoinMemcpyN(columnUpper_, numberColumns_, saveUpper);
4973
4974	int iRow, iColumn;
4975	// If wanted compute a reasonable dualBound_
4976	if (factor == COIN_DBL_MAX) {
4977	factor = 0.0;
4978	if (dualBound_ == 1.0e10) {
4979	// get largest scaled away from bound
4980	double largest = 1.0e-12;
4981	double largestScaled = 1.0e-12;
4982	int iRow;
4983	for (iRow = 0; iRow < numberRows_; iRow++) {
4984	double value = rowActivity_[iRow];
4985	double above = value - rowLower_[iRow];
4986	double below = rowUpper_[iRow] - value;
4987	if (above < 1.0e12) {
4988	largest = CoinMax(largest, above);
4989	}
4990	if (below < 1.0e12) {
4991	largest = CoinMax(largest, below);
4992	}
4993	if (rowScale_) {
4994	double multiplier = rowScale_[iRow];
4995	above *= multiplier;
4996	below *= multiplier;
4997	}
4998	if (above < 1.0e12) {
4999	largestScaled = CoinMax(largestScaled, above);
5000	}
5001	if (below < 1.0e12) {
5002	largestScaled = CoinMax(largestScaled, below);
5003	}
5004	}
5005
5006	int iColumn;
5007	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5008	double value = columnActivity_[iColumn];
5009	double above = value - columnLower_[iColumn];
5010	double below = columnUpper_[iColumn] - value;
5011	if (above < 1.0e12) {
5012	largest = CoinMax(largest, above);
5013	}
5014	if (below < 1.0e12) {
5015	largest = CoinMax(largest, below);
5016	}
5017	if (columnScale_) {
5018	double multiplier = 1.0 / columnScale_[iColumn];
5019	above *= multiplier;
5020	below *= multiplier;
5021	}
5022	if (above < 1.0e12) {
5023	largestScaled = CoinMax(largestScaled, above);
5024	}
5025	if (below < 1.0e12) {
5026	largestScaled = CoinMax(largestScaled, below);
5027	}
5028	}
5029	std::cout << "Largest (scaled) away from bound " << largestScaled
5030	<< " unscaled " << largest << std::endl;
5031	dualBound_ = CoinMax(1.0001e7, CoinMin(100.0 * largest, 1.00001e10));
5032	}
5033	}
5034
5035	// If wanted - tighten column bounds using solution
5036	if (factor) {
5037	double largest = 0.0;
5038	if (factor > 0.0) {
5039	assert (factor > 1.0);
5040	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5041	if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
5042	largest = CoinMax(largest, fabs(columnActivity_[iColumn]));
5043	}
5044	}
5045	largest *= factor;
5046	} else {
5047	// absolute
5048	largest = - factor;
5049	}
5050	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
5051	if (columnUpper_[iColumn] - columnLower_[iColumn] > tolerance) {
5052	columnUpper_[iColumn] = CoinMin(columnUpper_[iColumn], largest);
5053	columnLower_[iColumn] = CoinMax(columnLower_[iColumn], -largest);
5054	}
5055	}
5056	}
5057	#define MAXPASS 10
5058
5059	// Loop round seeing if we can tighten bounds
5060	// Would be faster to have a stack of possible rows
5061	// and we put altered rows back on stack
5062	int numberCheck = -1;
5063	while(numberChanged > numberCheck) {
5064
5065	numberChanged = 0; // Bounds tightened this pass
5066
5067	if (iPass == MAXPASS) break;
5068	iPass++;
5069
5070	for (iRow = 0; iRow < numberRows_; iRow++) {
5071
5072	if (rowLower_[iRow] > -large \|\| rowUpper_[iRow] < large) {
5073
5074	// possible row
5075	int infiniteUpper = 0;
5076	int infiniteLower = 0;
5077	double maximumUp = 0.0;
5078	double maximumDown = 0.0;
5079	double newBound;
5080	CoinBigIndex rStart = rowStart[iRow];
5081	CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
5082	CoinBigIndex j;
5083	// Compute possible lower and upper ranges
5084
5085	for (j = rStart; j < rEnd; ++j) {
5086	double value = element[j];
5087	iColumn = column[j];
5088	if (value > 0.0) {
5089	if (columnUpper_[iColumn] >= large) {
5090	++infiniteUpper;
5091	} else {
5092	maximumUp += columnUpper_[iColumn] * value;
5093	}
5094	if (columnLower_[iColumn] <= -large) {
5095	++infiniteLower;
5096	} else {
5097	maximumDown += columnLower_[iColumn] * value;
5098	}
5099	} else if (value < 0.0) {
5100	if (columnUpper_[iColumn] >= large) {
5101	++infiniteLower;
5102	} else {
5103	maximumDown += columnUpper_[iColumn] * value;
5104	}
5105	if (columnLower_[iColumn] <= -large) {
5106	++infiniteUpper;
5107	} else {
5108	maximumUp += columnLower_[iColumn] * value;
5109	}
5110	}
5111	}
5112	// Build in a margin of error
5113	maximumUp += 1.0e-8 * fabs(maximumUp);
5114	maximumDown -= 1.0e-8 * fabs(maximumDown);
5115	double maxUp = maximumUp + infiniteUpper * 1.0e31;
5116	double maxDown = maximumDown - infiniteLower * 1.0e31;
5117	if (maxUp <= rowUpper_[iRow] + tolerance &&
5118	maxDown >= rowLower_[iRow] - tolerance) {
5119
5120	// Row is redundant - make totally free
5121	// NO - can't do this for postsolve
5122	// rowLower_[iRow]=-COIN_DBL_MAX;
5123	// rowUpper_[iRow]=COIN_DBL_MAX;
5124	//printf("Redundant row in presolveX %d\n",iRow);
5125
5126	} else {
5127	if (maxUp < rowLower_[iRow] - 100.0 * tolerance \|\|
5128	maxDown > rowUpper_[iRow] + 100.0 * tolerance) {
5129	// problem is infeasible - exit at once
5130	numberInfeasible++;
5131	break;
5132	}
5133	double lower = rowLower_[iRow];
5134	double upper = rowUpper_[iRow];
5135	for (j = rStart; j < rEnd; ++j) {
5136	double value = element[j];
5137	iColumn = column[j];
5138	double nowLower = columnLower_[iColumn];
5139	double nowUpper = columnUpper_[iColumn];
5140	if (value > 0.0) {
5141	// positive value
5142	if (lower > -large) {
5143	if (!infiniteUpper) {
5144	assert(nowUpper < large2);
5145	newBound = nowUpper +
5146	(lower - maximumUp) / value;
5147	// relax if original was large
5148	if (fabs(maximumUp) > 1.0e8)
5149	newBound -= 1.0e-12 * fabs(maximumUp);
5150	} else if (infiniteUpper == 1 && nowUpper > large) {
5151	newBound = (lower - maximumUp) / value;
5152	// relax if original was large
5153	if (fabs(maximumUp) > 1.0e8)
5154	newBound -= 1.0e-