Context Navigation

ClpSimplex.cpp @ 1441

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