Context Navigation

source: trunk/ClpSimplex.cpp @ 372

Last change on this file since 372 was 372, checked in by forrest, 16 years ago
Moving scaling to model plus robustness fixes in simplex still messing with interior point
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 176.3 KB

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