Context Navigation

ClpSimplex.cpp @ 2393

Last change on this file since 2393 was 2393, checked in by forrest, 2 years ago
change problem status of 4 to 1 for bound conflict
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 436.0 KB

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