Context Navigation

source: branches/pre/ClpSimplex.cpp @ 182

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: