Context Navigation

ClpSimplex.cpp @ 1404

Last change on this file since 1404 was 1404, checked in by tkr, 12 years ago
Merging r1375:1393 from stable/BSP to stable/1.10
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 341.9 KB

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