Context Navigation

ClpSimplex.cpp @ 1785

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

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