Context Navigation

ClpSimplex.cpp @ 223

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

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

Context Navigation

source: branches/pre/ClpSimplex.cpp @ 223

Download in other formats: