Context Navigation

ClpSimplexOther.cpp @ 1784

Last change on this file since 1784 was 1784, checked in by forrest, 9 years ago
minor change for null objective
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 259.2 KB

Line
1	/* $Id: ClpSimplexOther.cpp 1784 2011-08-22 08:32:06Z forrest $ */
2	// Copyright (C) 2004, International Business Machines
3	// Corporation and others. All Rights Reserved.
4	// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6	#include "CoinPragma.hpp"
7
8	#include <math.h>
9
10	#include "CoinHelperFunctions.hpp"
11	#include "ClpSimplexOther.hpp"
12	#include "ClpSimplexDual.hpp"
13	#include "ClpSimplexPrimal.hpp"
14	#include "ClpEventHandler.hpp"
15	#include "ClpHelperFunctions.hpp"
16	#include "ClpFactorization.hpp"
17	#include "ClpDualRowDantzig.hpp"
18	#include "ClpNonLinearCost.hpp"
19	#include "ClpDynamicMatrix.hpp"
20	#include "CoinPackedMatrix.hpp"
21	#include "CoinIndexedVector.hpp"
22	#include "CoinBuild.hpp"
23	#include "CoinMpsIO.hpp"
24	#include "CoinFloatEqual.hpp"
25	#include "ClpMessage.hpp"
26	#include <cfloat>
27	#include <cassert>
28	#include <string>
29	#include <stdio.h>
30	#include <iostream>
31	/* Dual ranging.
32	This computes increase/decrease in cost for each given variable and corresponding
33	sequence numbers which would change basis. Sequence numbers are 0..numberColumns
34	and numberColumns.. for artificials/slacks.
35	For non-basic variables the sequence number will be that of the non-basic variables.
36
37	Up to user to provide correct length arrays.
38
39	*/
40	void ClpSimplexOther::dualRanging(int numberCheck, const int * which,
41	double * costIncreased, int * sequenceIncreased,
42	double * costDecreased, int * sequenceDecreased,
43	double * valueIncrease, double * valueDecrease)
44	{
45	rowArray_[1]->clear();
46	columnArray_[1]->clear();
47	// long enough for rows+columns
48	assert(rowArray_[3]->capacity() >= numberRows_ + numberColumns_);
49	rowArray_[3]->clear();
50	int * backPivot = rowArray_[3]->getIndices();
51	int i;
52	for ( i = 0; i < numberRows_ + numberColumns_; i++) {
53	backPivot[i] = -1;
54	}
55	for (i = 0; i < numberRows_; i++) {
56	int iSequence = pivotVariable_[i];
57	backPivot[iSequence] = i;
58	}
59	// dualTolerance may be zero if from CBC. In fact use that fact
60	bool inCBC = !dualTolerance_;
61	if (inCBC)
62	assert (integerType_);
63	dualTolerance_ = dblParam_[ClpDualTolerance];
64	double * arrayX = rowArray_[0]->denseVector();
65	for ( i = 0; i < numberCheck; i++) {
66	rowArray_[0]->clear();
67	//rowArray_[0]->checkClear();
68	//rowArray_[1]->checkClear();
69	//columnArray_[1]->checkClear();
70	columnArray_[0]->clear();
71	//columnArray_[0]->checkClear();
72	int iSequence = which[i];
73	if (iSequence < 0) {
74	costIncreased[i] = 0.0;
75	sequenceIncreased[i] = -1;
76	costDecreased[i] = 0.0;
77	sequenceDecreased[i] = -1;
78	continue;
79	}
80	double costIncrease = COIN_DBL_MAX;
81	double costDecrease = COIN_DBL_MAX;
82	int sequenceIncrease = -1;
83	int sequenceDecrease = -1;
84	if (valueIncrease) {
85	assert (valueDecrease);
86	valueIncrease[i] = iSequence < numberColumns_ ? columnActivity_[iSequence] : rowActivity_[iSequence-numberColumns_];
87	valueDecrease[i] = valueIncrease[i];
88	}
89
90	switch(getStatus(iSequence)) {
91
92	case basic: {
93	// non-trvial
94	// Get pivot row
95	int iRow = backPivot[iSequence];
96	assert (iRow >= 0);
97	double plusOne = 1.0;
98	rowArray_[0]->createPacked(1, &iRow, &plusOne);
99	factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
100	// put row of tableau in rowArray[0] and columnArray[0]
101	matrix_->transposeTimes(this, -1.0,
102	rowArray_[0], columnArray_[1], columnArray_[0]);
103	double alphaIncrease;
104	double alphaDecrease;
105	// do ratio test up and down
106	checkDualRatios(rowArray_[0], columnArray_[0], costIncrease, sequenceIncrease, alphaIncrease,
107	costDecrease, sequenceDecrease, alphaDecrease);
108	if (!inCBC) {
109	if (valueIncrease) {
110	if (sequenceIncrease >= 0)
111	valueIncrease[i] = primalRanging1(sequenceIncrease, iSequence);
112	if (sequenceDecrease >= 0)
113	valueDecrease[i] = primalRanging1(sequenceDecrease, iSequence);
114	}
115	} else {
116	int number = rowArray_[0]->getNumElements();
117	double scale2 = 0.0;
118	int j;
119	for (j = 0; j < number; j++) {
120	scale2 += arrayX[j] * arrayX[j];
121	}
122	scale2 = 1.0 / sqrt(scale2);
123	//valueIncrease[i] = scale2;
124	if (sequenceIncrease >= 0) {
125	double djValue = dj_[sequenceIncrease];
126	if (fabs(djValue) > 10.0 * dualTolerance_) {
127	// we are going to use for cutoff so be exact
128	costIncrease = fabs(djValue / alphaIncrease);
129	/* Not sure this is good idea as I don't think correct e.g.
130	suppose a continuous variable has dj slightly greater. */
131	if(false && sequenceIncrease < numberColumns_ && integerType_[sequenceIncrease]) {
132	// can improve
133	double movement = (columnScale_ == NULL) ? 1.0 :
134	rhsScale_ * inverseColumnScale_[sequenceIncrease];
135	costIncrease = CoinMax(fabs(djValue * movement), costIncrease);
136	}
137	} else {
138	costIncrease = 0.0;
139	}
140	}
141	if (sequenceDecrease >= 0) {
142	double djValue = dj_[sequenceDecrease];
143	if (fabs(djValue) > 10.0 * dualTolerance_) {
144	// we are going to use for cutoff so be exact
145	costDecrease = fabs(djValue / alphaDecrease);
146	if(sequenceDecrease < numberColumns_ && integerType_[sequenceDecrease]) {
147	// can improve
148	double movement = (columnScale_ == NULL) ? 1.0 :
149	rhsScale_ * inverseColumnScale_[sequenceDecrease];
150	costDecrease = CoinMax(fabs(djValue * movement), costDecrease);
151	}
152	} else {
153	costDecrease = 0.0;
154	}
155	}
156	costIncrease *= scale2;
157	costDecrease *= scale2;
158	}
159	}
160	break;
161	case isFixed:
162	break;
163	case isFree:
164	case superBasic:
165	costIncrease = 0.0;
166	costDecrease = 0.0;
167	sequenceIncrease = iSequence;
168	sequenceDecrease = iSequence;
169	break;
170	case atUpperBound:
171	costIncrease = CoinMax(0.0, -dj_[iSequence]);
172	sequenceIncrease = iSequence;
173	if (valueIncrease)
174	valueIncrease[i] = primalRanging1(iSequence, iSequence);
175	break;
176	case atLowerBound:
177	costDecrease = CoinMax(0.0, dj_[iSequence]);
178	sequenceDecrease = iSequence;
179	if (valueIncrease)
180	valueDecrease[i] = primalRanging1(iSequence, iSequence);
181	break;
182	}
183	double scaleFactor;
184	if (rowScale_) {
185	if (iSequence < numberColumns_)
186	scaleFactor = 1.0 / (objectiveScale_ * columnScale_[iSequence]);
187	else
188	scaleFactor = rowScale_[iSequence-numberColumns_] / objectiveScale_;
189	} else {
190	scaleFactor = 1.0 / objectiveScale_;
191	}
192	if (costIncrease < 1.0e30)
193	costIncrease *= scaleFactor;
194	if (costDecrease < 1.0e30)
195	costDecrease *= scaleFactor;
196	if (optimizationDirection_ == 1.0) {
197	costIncreased[i] = costIncrease;
198	sequenceIncreased[i] = sequenceIncrease;
199	costDecreased[i] = costDecrease;
200	sequenceDecreased[i] = sequenceDecrease;
201	} else if (optimizationDirection_ == -1.0) {
202	costIncreased[i] = costDecrease;
203	sequenceIncreased[i] = sequenceDecrease;
204	costDecreased[i] = costIncrease;
205	sequenceDecreased[i] = sequenceIncrease;
206	if (valueIncrease) {
207	double temp = valueIncrease[i];
208	valueIncrease[i] = valueDecrease[i];
209	valueDecrease[i] = temp;
210	}
211	} else if (optimizationDirection_ == 0.0) {
212	// !!!!!! ???
213	costIncreased[i] = COIN_DBL_MAX;
214	sequenceIncreased[i] = -1;
215	costDecreased[i] = COIN_DBL_MAX;
216	sequenceDecreased[i] = -1;
217	} else {
218	abort();
219	}
220	}
221	rowArray_[0]->clear();
222	//rowArray_[1]->clear();
223	//columnArray_[1]->clear();
224	columnArray_[0]->clear();
225	//rowArray_[3]->clear();
226	if (!optimizationDirection_)
227	printf("*** ????? Ranging with zero optimization costs\n");
228	}
229	/*
230	Row array has row part of pivot row
231	Column array has column part.
232	This is used in dual ranging
233	*/
234	void
235	ClpSimplexOther::checkDualRatios(CoinIndexedVector * rowArray,
236	CoinIndexedVector * columnArray,
237	double & costIncrease, int & sequenceIncrease, double & alphaIncrease,
238	double & costDecrease, int & sequenceDecrease, double & alphaDecrease)
239	{
240	double acceptablePivot = 1.0e-9;
241	double * work;
242	int number;
243	int * which;
244	int iSection;
245
246	double thetaDown = 1.0e31;
247	double thetaUp = 1.0e31;
248	int sequenceDown = -1;
249	int sequenceUp = -1;
250	double alphaDown = 0.0;
251	double alphaUp = 0.0;
252
253	int addSequence;
254
255	for (iSection = 0; iSection < 2; iSection++) {
256
257	int i;
258	if (!iSection) {
259	work = rowArray->denseVector();
260	number = rowArray->getNumElements();
261	which = rowArray->getIndices();
262	addSequence = numberColumns_;
263	} else {
264	work = columnArray->denseVector();
265	number = columnArray->getNumElements();
266	which = columnArray->getIndices();
267	addSequence = 0;
268	}
269
270	for (i = 0; i < number; i++) {
271	int iSequence = which[i];
272	int iSequence2 = iSequence + addSequence;
273	double alpha = work[i];
274	if (fabs(alpha) < acceptablePivot)
275	continue;
276	double oldValue = dj_[iSequence2];
277
278	switch(getStatus(iSequence2)) {
279
280	case basic:
281	break;
282	case ClpSimplex::isFixed:
283	break;
284	case isFree:
285	case superBasic:
286	// treat dj as if zero
287	thetaDown = 0.0;
288	thetaUp = 0.0;
289	sequenceDown = iSequence2;
290	sequenceUp = iSequence2;
291	break;
292	case atUpperBound:
293	if (alpha > 0.0) {
294	// test up
295	if (oldValue + thetaUp * alpha > dualTolerance_) {
296	thetaUp = (dualTolerance_ - oldValue) / alpha;
297	sequenceUp = iSequence2;
298	alphaUp = alpha;
299	}
300	} else {
301	// test down
302	if (oldValue - thetaDown * alpha > dualTolerance_) {
303	thetaDown = -(dualTolerance_ - oldValue) / alpha;
304	sequenceDown = iSequence2;
305	alphaDown = alpha;
306	}
307	}
308	break;
309	case atLowerBound:
310	if (alpha < 0.0) {
311	// test up
312	if (oldValue + thetaUp * alpha < - dualTolerance_) {
313	thetaUp = -(dualTolerance_ + oldValue) / alpha;
314	sequenceUp = iSequence2;
315	alphaUp = alpha;
316	}
317	} else {
318	// test down
319	if (oldValue - thetaDown * alpha < -dualTolerance_) {
320	thetaDown = (dualTolerance_ + oldValue) / alpha;
321	sequenceDown = iSequence2;
322	alphaDown = alpha;
323	}
324	}
325	break;
326	}
327	}
328	}
329	if (sequenceUp >= 0) {
330	costIncrease = thetaUp;
331	sequenceIncrease = sequenceUp;
332	alphaIncrease = alphaUp;
333	}
334	if (sequenceDown >= 0) {
335	costDecrease = thetaDown;
336	sequenceDecrease = sequenceDown;
337	alphaDecrease = alphaDown;
338	}
339	}
340	/** Primal ranging.
341	This computes increase/decrease in value for each given variable and corresponding
342	sequence numbers which would change basis. Sequence numbers are 0..numberColumns
343	and numberColumns.. for artificials/slacks.
344	For basic variables the sequence number will be that of the basic variables.
345
346	Up to user to provide correct length arrays.
347
348	When here - guaranteed optimal
349	*/
350	void
351	ClpSimplexOther::primalRanging(int numberCheck, const int * which,
352	double * valueIncreased, int * sequenceIncreased,
353	double * valueDecreased, int * sequenceDecreased)
354	{
355	rowArray_[0]->clear();
356	rowArray_[1]->clear();
357	lowerIn_ = -COIN_DBL_MAX;
358	upperIn_ = COIN_DBL_MAX;
359	valueIn_ = 0.0;
360	for ( int i = 0; i < numberCheck; i++) {
361	int iSequence = which[i];
362	double valueIncrease = COIN_DBL_MAX;
363	double valueDecrease = COIN_DBL_MAX;
364	int sequenceIncrease = -1;
365	int sequenceDecrease = -1;
366
367	switch(getStatus(iSequence)) {
368
369	case basic:
370	case isFree:
371	case superBasic:
372	// Easy
373	valueDecrease = CoinMax(0.0, upper_[iSequence] - solution_[iSequence]);
374	valueIncrease = CoinMax(0.0, solution_[iSequence] - lower_[iSequence]);
375	sequenceDecrease = iSequence;
376	sequenceIncrease = iSequence;
377	break;
378	case isFixed:
379	case atUpperBound:
380	case atLowerBound: {
381	// Non trivial
382	// Other bound is ignored
383	unpackPacked(rowArray_[1], iSequence);
384	factorization_->updateColumn(rowArray_[2], rowArray_[1]);
385	// Get extra rows
386	matrix_->extendUpdated(this, rowArray_[1], 0);
387	// do ratio test
388	checkPrimalRatios(rowArray_[1], 1);
389	if (pivotRow_ >= 0) {
390	valueIncrease = theta_;
391	sequenceIncrease = pivotVariable_[pivotRow_];
392	}
393	checkPrimalRatios(rowArray_[1], -1);
394	if (pivotRow_ >= 0) {
395	valueDecrease = theta_;
396	sequenceDecrease = pivotVariable_[pivotRow_];
397	}
398	rowArray_[1]->clear();
399	}
400	break;
401	}
402	double scaleFactor;
403	if (rowScale_) {
404	if (iSequence < numberColumns_)
405	scaleFactor = columnScale_[iSequence] / rhsScale_;
406	else
407	scaleFactor = 1.0 / (rowScale_[iSequence-numberColumns_] * rhsScale_);
408	} else {
409	scaleFactor = 1.0 / rhsScale_;
410	}
411	if (valueIncrease < 1.0e30)
412	valueIncrease *= scaleFactor;
413	else
414	valueIncrease = COIN_DBL_MAX;
415	if (valueDecrease < 1.0e30)
416	valueDecrease *= scaleFactor;
417	else
418	valueDecrease = COIN_DBL_MAX;
419	valueIncreased[i] = valueIncrease;
420	sequenceIncreased[i] = sequenceIncrease;
421	valueDecreased[i] = valueDecrease;
422	sequenceDecreased[i] = sequenceDecrease;
423	}
424	}
425	// Returns new value of whichOther when whichIn enters basis
426	double
427	ClpSimplexOther::primalRanging1(int whichIn, int whichOther)
428	{
429	rowArray_[0]->clear();
430	rowArray_[1]->clear();
431	int iSequence = whichIn;
432	double newValue = solution_[whichOther];
433	double alphaOther = 0.0;
434	Status status = getStatus(iSequence);
435	assert (status == atLowerBound \|\| status == atUpperBound);
436	int wayIn = (status == atLowerBound) ? 1 : -1;
437
438	switch(getStatus(iSequence)) {
439
440	case basic:
441	case isFree:
442	case superBasic:
443	assert (whichIn == whichOther);
444	// Easy
445	newValue = wayIn > 0 ? upper_[iSequence] : lower_[iSequence];
446	break;
447	case isFixed:
448	case atUpperBound:
449	case atLowerBound:
450	// Non trivial
451	{
452	// Other bound is ignored
453	unpackPacked(rowArray_[1], iSequence);
454	factorization_->updateColumn(rowArray_[2], rowArray_[1]);
455	// Get extra rows
456	matrix_->extendUpdated(this, rowArray_[1], 0);
457	// do ratio test
458	double acceptablePivot = 1.0e-7;
459	double * work = rowArray_[1]->denseVector();
460	int number = rowArray_[1]->getNumElements();
461	int * which = rowArray_[1]->getIndices();
462
463	// we may need to swap sign
464	double way = wayIn;
465	double theta = 1.0e30;
466	for (int iIndex = 0; iIndex < number; iIndex++) {
467
468	int iRow = which[iIndex];
469	double alpha = work[iIndex] * way;
470	int iPivot = pivotVariable_[iRow];
471	if (iPivot == whichOther) {
472	alphaOther = alpha;
473	continue;
474	}
475	double oldValue = solution_[iPivot];
476	if (fabs(alpha) > acceptablePivot) {
477	if (alpha > 0.0) {
478	// basic variable going towards lower bound
479	double bound = lower_[iPivot];
480	oldValue -= bound;
481	if (oldValue - theta * alpha < 0.0) {
482	theta = CoinMax(0.0, oldValue / alpha);
483	}
484	} else {
485	// basic variable going towards upper bound
486	double bound = upper_[iPivot];
487	oldValue = oldValue - bound;
488	if (oldValue - theta * alpha > 0.0) {
489	theta = CoinMax(0.0, oldValue / alpha);
490	}
491	}
492	}
493	}
494	if (whichIn != whichOther) {
495	if (theta < 1.0e30)
496	newValue -= theta * alphaOther;
497	else
498	newValue = alphaOther > 0.0 ? -1.0e30 : 1.0e30;
499	} else {
500	newValue += theta * wayIn;
501	}
502	}
503	rowArray_[1]->clear();
504	break;
505	}
506	double scaleFactor;
507	if (rowScale_) {
508	if (whichOther < numberColumns_)
509	scaleFactor = columnScale_[whichOther] / rhsScale_;
510	else
511	scaleFactor = 1.0 / (rowScale_[whichOther-numberColumns_] * rhsScale_);
512	} else {
513	scaleFactor = 1.0 / rhsScale_;
514	}
515	if (newValue < 1.0e29)
516	if (newValue > -1.0e29)
517	newValue *= scaleFactor;
518	else
519	newValue = -COIN_DBL_MAX;
520	else
521	newValue = COIN_DBL_MAX;
522	return newValue;
523	}
524	/*
525	Row array has pivot column
526	This is used in primal ranging
527	*/
528	void
529	ClpSimplexOther::checkPrimalRatios(CoinIndexedVector * rowArray,
530	int direction)
531	{
532	// sequence stays as row number until end
533	pivotRow_ = -1;
534	double acceptablePivot = 1.0e-7;
535	double * work = rowArray->denseVector();
536	int number = rowArray->getNumElements();
537	int * which = rowArray->getIndices();
538
539	// we need to swap sign if going down
540	double way = direction;
541	theta_ = 1.0e30;
542	for (int iIndex = 0; iIndex < number; iIndex++) {
543
544	int iRow = which[iIndex];
545	double alpha = work[iIndex] * way;
546	int iPivot = pivotVariable_[iRow];
547	double oldValue = solution_[iPivot];
548	if (fabs(alpha) > acceptablePivot) {
549	if (alpha > 0.0) {
550	// basic variable going towards lower bound
551	double bound = lower_[iPivot];
552	oldValue -= bound;
553	if (oldValue - theta_ * alpha < 0.0) {
554	pivotRow_ = iRow;
555	theta_ = CoinMax(0.0, oldValue / alpha);
556	}
557	} else {
558	// basic variable going towards upper bound
559	double bound = upper_[iPivot];
560	oldValue = oldValue - bound;
561	if (oldValue - theta_ * alpha > 0.0) {
562	pivotRow_ = iRow;
563	theta_ = CoinMax(0.0, oldValue / alpha);
564	}
565	}
566	}
567	}
568	}
569	/* Write the basis in MPS format to the specified file.
570	If writeValues true writes values of structurals
571	(and adds VALUES to end of NAME card)
572
573	Row and column names may be null.
574	formatType is
575	<ul>
576	<li> 0 - normal
577	<li> 1 - extra accuracy
578	<li> 2 - IEEE hex (later)
579	</ul>
580
581	Returns non-zero on I/O error
582
583	This is based on code contributed by Thorsten Koch
584	*/
585	int
586	ClpSimplexOther::writeBasis(const char *filename,
587	bool writeValues,
588	int formatType) const
589	{
590	formatType = CoinMax(0, formatType);
591	formatType = CoinMin(2, formatType);
592	if (!writeValues)
593	formatType = 0;
594	// See if INTEL if IEEE
595	if (formatType == 2) {
596	// test intel here and add 1 if not intel
597	double value = 1.0;
598	char x[8];
599	memcpy(x, &value, 8);
600	if (x[0] == 63) {
601	formatType ++; // not intel
602	} else {
603	assert (x[0] == 0);
604	}
605	}
606
607	char number[20];
608	FILE * fp = fopen(filename, "w");
609	if (!fp)
610	return -1;
611
612	// NAME card
613
614	if (strcmp(strParam_[ClpProbName].c_str(), "") == 0) {
615	fprintf(fp, "NAME BLANK ");
616	} else {
617	fprintf(fp, "NAME %s ", strParam_[ClpProbName].c_str());
618	}
619	if (formatType >= 2)
620	fprintf(fp, "FREEIEEE");
621	else if (writeValues)
622	fprintf(fp, "VALUES");
623	// finish off name
624	fprintf(fp, "\n");
625	int iRow = 0;
626	for(int iColumn = 0; iColumn < numberColumns_; iColumn++) {
627	bool printit = false;
628	if( getColumnStatus(iColumn) == ClpSimplex::basic) {
629	printit = true;
630	// Find non basic row
631	for(; iRow < numberRows_; iRow++) {
632	if (getRowStatus(iRow) != ClpSimplex::basic)
633	break;
634	}
635	if (lengthNames_) {
636	if (iRow != numberRows_) {
637	fprintf(fp, " %s %-8s %s",
638	getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
639	columnNames_[iColumn].c_str(),
640	rowNames_[iRow].c_str());
641	iRow++;
642	} else {
643	// Allow for too many basics!
644	fprintf(fp, " BS %-8s ",
645	columnNames_[iColumn].c_str());
646	// Dummy row name if values
647	if (writeValues)
648	fprintf(fp, " _dummy_");
649	}
650	} else {
651	// no names
652	if (iRow != numberRows_) {
653	fprintf(fp, " %s C%7.7d R%7.7d",
654	getRowStatus(iRow) == ClpSimplex::atUpperBound ? "XU" : "XL",
655	iColumn, iRow);
656	iRow++;
657	} else {
658	// Allow for too many basics!
659	fprintf(fp, " BS C%7.7d", iColumn);
660	// Dummy row name if values
661	if (writeValues)
662	fprintf(fp, " _dummy_");
663	}
664	}
665	} else {
666	if( getColumnStatus(iColumn) == ClpSimplex::atUpperBound) {
667	printit = true;
668	if (lengthNames_)
669	fprintf(fp, " UL %s", columnNames_[iColumn].c_str());
670	else
671	fprintf(fp, " UL C%7.7d", iColumn);
672	// Dummy row name if values
673	if (writeValues)
674	fprintf(fp, " _dummy_");
675	}
676	}
677	if (printit && writeValues) {
678	// add value
679	CoinConvertDouble(0, formatType, columnActivity_[iColumn], number);
680	fprintf(fp, " %s", number);
681	}
682	if (printit)
683	fprintf(fp, "\n");
684	}
685	fprintf(fp, "ENDATA\n");
686	fclose(fp);
687	return 0;
688	}
689	// Read a basis from the given filename
690	int
691	ClpSimplexOther::readBasis(const char *fileName)
692	{
693	int status = 0;
694	bool canOpen = false;
695	if (!strcmp(fileName, "-") \|\| !strcmp(fileName, "stdin")) {
696	// stdin
697	canOpen = true;
698	} else {
699	FILE *fp = fopen(fileName, "r");
700	if (fp) {
701	// can open - lets go for it
702	fclose(fp);
703	canOpen = true;
704	} else {
705	handler_->message(CLP_UNABLE_OPEN, messages_)
706	<< fileName << CoinMessageEol;
707	return -1;
708	}
709	}
710	CoinMpsIO m;
711	m.passInMessageHandler(handler_);
712	*m.messagesPointer() = coinMessages();
713	bool savePrefix = m.messageHandler()->prefix();
714	m.messageHandler()->setPrefix(handler_->prefix());
715	status = m.readBasis(fileName, "", columnActivity_, status_ + numberColumns_,
716	status_,
717	columnNames_, numberColumns_,
718	rowNames_, numberRows_);
719	m.messageHandler()->setPrefix(savePrefix);
720	if (status >= 0) {
721	if (!status) {
722	// set values
723	int iColumn, iRow;
724	for (iRow = 0; iRow < numberRows_; iRow++) {
725	if (getRowStatus(iRow) == atLowerBound)
726	rowActivity_[iRow] = rowLower_[iRow];
727	else if (getRowStatus(iRow) == atUpperBound)
728	rowActivity_[iRow] = rowUpper_[iRow];
729	}
730	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
731	if (getColumnStatus(iColumn) == atLowerBound)
732	columnActivity_[iColumn] = columnLower_[iColumn];
733	else if (getColumnStatus(iColumn) == atUpperBound)
734	columnActivity_[iColumn] = columnUpper_[iColumn];
735	}
736	} else {
737	memset(rowActivity_, 0, numberRows_ * sizeof(double));
738	matrix_->times(-1.0, columnActivity_, rowActivity_);
739	}
740	} else {
741	// errors
742	handler_->message(CLP_IMPORT_ERRORS, messages_)
743	<< status << fileName << CoinMessageEol;
744	}
745	return status;
746	}
747	/* Creates dual of a problem if looks plausible
748	(defaults will always create model)
749	fractionRowRanges is fraction of rows allowed to have ranges
750	fractionColumnRanges is fraction of columns allowed to have ranges
751	*/
752	ClpSimplex *
753	ClpSimplexOther::dualOfModel(double fractionRowRanges, double fractionColumnRanges) const
754	{
755	const ClpSimplex * model2 = static_cast<const ClpSimplex *> (this);
756	bool changed = false;
757	int numberChanged = 0;
758	int iColumn;
759	// check if we need to change bounds to rows
760	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
761	if (columnUpper_[iColumn] < 1.0e20 &&
762	columnLower_[iColumn] > -1.0e20) {
763	changed = true;
764	numberChanged++;
765	}
766	}
767	int iRow;
768	int numberExtraRows = 0;
769	if (numberChanged <= fractionColumnRanges * numberColumns_) {
770	for (iRow = 0; iRow < numberRows_; iRow++) {
771	if (rowLower_[iRow] > -1.0e20 &&
772	rowUpper_[iRow] < 1.0e20) {
773	if (rowUpper_[iRow] != rowLower_[iRow])
774	numberExtraRows++;
775	}
776	}
777	if (numberExtraRows > fractionRowRanges * numberRows_)
778	return NULL;
779	} else {
780	return NULL;
781	}
782	if (changed) {
783	ClpSimplex * model3 = new ClpSimplex(*model2);
784	CoinBuild build;
785	double one = 1.0;
786	int numberColumns = model3->numberColumns();
787	const double * columnLower = model3->columnLower();
788	const double * columnUpper = model3->columnUpper();
789	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
790	if (columnUpper[iColumn] < 1.0e20 &&
791	columnLower[iColumn] > -1.0e20) {
792	if (fabs(columnLower[iColumn]) < fabs(columnUpper[iColumn])) {
793	double value = columnUpper[iColumn];
794	model3->setColumnUpper(iColumn, COIN_DBL_MAX);
795	build.addRow(1, &iColumn, &one, -COIN_DBL_MAX, value);
796	} else {
797	double value = columnLower[iColumn];
798	model3->setColumnLower(iColumn, -COIN_DBL_MAX);
799	build.addRow(1, &iColumn, &one, value, COIN_DBL_MAX);
800	}
801	}
802	}
803	model3->addRows(build);
804	model2 = model3;
805	}
806	int numberColumns = model2->numberColumns();
807	const double * columnLower = model2->columnLower();
808	const double * columnUpper = model2->columnUpper();
809	int numberRows = model2->numberRows();
810	double * rowLower = CoinCopyOfArray(model2->rowLower(), numberRows);
811	double * rowUpper = CoinCopyOfArray(model2->rowUpper(), numberRows);
812
813	const double * objective = model2->objective();
814	CoinPackedMatrix * matrix = model2->matrix();
815	// get transpose
816	CoinPackedMatrix rowCopy = *matrix;
817	const int * row = matrix->getIndices();
818	const int * columnLength = matrix->getVectorLengths();
819	const CoinBigIndex * columnStart = matrix->getVectorStarts();
820	const double * elementByColumn = matrix->getElements();
821	double objOffset = 0.0;
822	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
823	double offset = 0.0;
824	double objValue = optimizationDirection_ * objective[iColumn];
825	if (columnUpper[iColumn] > 1.0e20) {
826	if (columnLower[iColumn] > -1.0e20)
827	offset = columnLower[iColumn];
828	} else if (columnLower[iColumn] < -1.0e20) {
829	offset = columnUpper[iColumn];
830	} else {
831	// taken care of before
832	abort();
833	}
834	if (offset) {
835	objOffset += offset * objValue;
836	for (CoinBigIndex j = columnStart[iColumn];
837	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
838	int iRow = row[j];
839	if (rowLower[iRow] > -1.0e20)
840	rowLower[iRow] -= offset * elementByColumn[j];
841	if (rowUpper[iRow] < 1.0e20)
842	rowUpper[iRow] -= offset * elementByColumn[j];
843	}
844	}
845	}
846	int * which = new int[numberRows+numberExtraRows];
847	rowCopy.reverseOrdering();
848	rowCopy.transpose();
849	double * fromRowsLower = new double[numberRows+numberExtraRows];
850	double * fromRowsUpper = new double[numberRows+numberExtraRows];
851	double * newObjective = new double[numberRows+numberExtraRows];
852	double * fromColumnsLower = new double[numberColumns];
853	double * fromColumnsUpper = new double[numberColumns];
854	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
855	double objValue = optimizationDirection_ * objective[iColumn];
856	// Offset is already in
857	if (columnUpper[iColumn] > 1.0e20) {
858	if (columnLower[iColumn] > -1.0e20) {
859	fromColumnsLower[iColumn] = -COIN_DBL_MAX;
860	fromColumnsUpper[iColumn] = objValue;
861	} else {
862	// free
863	fromColumnsLower[iColumn] = objValue;
864	fromColumnsUpper[iColumn] = objValue;
865	}
866	} else if (columnLower[iColumn] < -1.0e20) {
867	fromColumnsLower[iColumn] = objValue;
868	fromColumnsUpper[iColumn] = COIN_DBL_MAX;
869	} else {
870	abort();
871	}
872	}
873	int kRow = 0;
874	int kExtraRow = numberRows;
875	for (iRow = 0; iRow < numberRows; iRow++) {
876	if (rowLower[iRow] < -1.0e20) {
877	assert (rowUpper[iRow] < 1.0e20);
878	newObjective[kRow] = -rowUpper[iRow];
879	fromRowsLower[kRow] = -COIN_DBL_MAX;
880	fromRowsUpper[kRow] = 0.0;
881	which[kRow] = iRow;
882	kRow++;
883	} else if (rowUpper[iRow] > 1.0e20) {
884	newObjective[kRow] = -rowLower[iRow];
885	fromRowsLower[kRow] = 0.0;
886	fromRowsUpper[kRow] = COIN_DBL_MAX;
887	which[kRow] = iRow;
888	kRow++;
889	} else {
890	if (rowUpper[iRow] == rowLower[iRow]) {
891	newObjective[kRow] = -rowLower[iRow];
892	fromRowsLower[kRow] = -COIN_DBL_MAX;;
893	fromRowsUpper[kRow] = COIN_DBL_MAX;
894	which[kRow] = iRow;
895	kRow++;
896	} else {
897	// range
898	newObjective[kRow] = -rowUpper[iRow];
899	fromRowsLower[kRow] = -COIN_DBL_MAX;
900	fromRowsUpper[kRow] = 0.0;
901	which[kRow] = iRow;
902	kRow++;
903	newObjective[kExtraRow] = -rowLower[iRow];
904	fromRowsLower[kExtraRow] = 0.0;
905	fromRowsUpper[kExtraRow] = COIN_DBL_MAX;
906	which[kExtraRow] = iRow;
907	kExtraRow++;
908	}
909	}
910	}
911	if (numberExtraRows) {
912	CoinPackedMatrix newCopy;
913	newCopy.setExtraGap(0.0);
914	newCopy.setExtraMajor(0.0);
915	newCopy.submatrixOfWithDuplicates(rowCopy, kExtraRow, which);
916	rowCopy = newCopy;
917	}
918	ClpSimplex * modelDual = new ClpSimplex();
919	modelDual->loadProblem(rowCopy, fromRowsLower, fromRowsUpper, newObjective,
920	fromColumnsLower, fromColumnsUpper);
921	modelDual->setObjectiveOffset(objOffset);
922	modelDual->setDualBound(model2->dualBound());
923	modelDual->setInfeasibilityCost(model2->infeasibilityCost());
924	modelDual->setDualTolerance(model2->dualTolerance());
925	modelDual->setPrimalTolerance(model2->primalTolerance());
926	modelDual->setPerturbation(model2->perturbation());
927	modelDual->setSpecialOptions(model2->specialOptions());
928	modelDual->setMoreSpecialOptions(model2->moreSpecialOptions());
929	delete [] fromRowsLower;
930	delete [] fromRowsUpper;
931	delete [] fromColumnsLower;
932	delete [] fromColumnsUpper;
933	delete [] newObjective;
934	delete [] which;
935	delete [] rowLower;
936	delete [] rowUpper;
937	if (changed)
938	delete model2;
939	modelDual->createStatus();
940	return modelDual;
941	}
942	// Restores solution from dualized problem
943	int
944	ClpSimplexOther::restoreFromDual(const ClpSimplex * dualProblem)
945	{
946	int returnCode = 0;;
947	createStatus();
948	// Number of rows in dual problem was original number of columns
949	assert (numberColumns_ == dualProblem->numberRows());
950	// If slack on d-row basic then column at bound otherwise column basic
951	// If d-column basic then rhs tight
952	int numberBasic = 0;
953	int iRow, iColumn = 0;
954	// Get number of extra rows from ranges
955	int numberExtraRows = 0;
956	for (iRow = 0; iRow < numberRows_; iRow++) {
957	if (rowLower_[iRow] > -1.0e20 &&
958	rowUpper_[iRow] < 1.0e20) {
959	if (rowUpper_[iRow] != rowLower_[iRow])
960	numberExtraRows++;
961	}
962	}
963	const double * objective = this->objective();
964	const double * dualDual = dualProblem->dualRowSolution();
965	const double * dualDj = dualProblem->dualColumnSolution();
966	const double * dualSol = dualProblem->primalColumnSolution();
967	const double * dualActs = dualProblem->primalRowSolution();
968	#if 0
969	ClpSimplex thisCopy = *this;
970	thisCopy.dual(); // for testing
971	const double * primalDual = thisCopy.dualRowSolution();
972	const double * primalDj = thisCopy.dualColumnSolution();
973	const double * primalSol = thisCopy.primalColumnSolution();
974	const double * primalActs = thisCopy.primalRowSolution();
975	char ss[] = {'F', 'B', 'U', 'L', 'S', 'F'};
976	printf ("Dual problem row info %d rows\n", dualProblem->numberRows());
977	for (iRow = 0; iRow < dualProblem->numberRows(); iRow++)
978	printf("%d at %c primal %g dual %g\n",
979	iRow, ss[dualProblem->getRowStatus(iRow)],
980	dualActs[iRow], dualDual[iRow]);
981	printf ("Dual problem column info %d columns\n", dualProblem->numberColumns());
982	for (iColumn = 0; iColumn < dualProblem->numberColumns(); iColumn++)
983	printf("%d at %c primal %g dual %g\n",
984	iColumn, ss[dualProblem->getColumnStatus(iColumn)],
985	dualSol[iColumn], dualDj[iColumn]);
986	printf ("Primal problem row info %d rows\n", thisCopy.numberRows());
987	for (iRow = 0; iRow < thisCopy.numberRows(); iRow++)
988	printf("%d at %c primal %g dual %g\n",
989	iRow, ss[thisCopy.getRowStatus(iRow)],
990	primalActs[iRow], primalDual[iRow]);
991	printf ("Primal problem column info %d columns\n", thisCopy.numberColumns());
992	for (iColumn = 0; iColumn < thisCopy.numberColumns(); iColumn++)
993	printf("%d at %c primal %g dual %g\n",
994	iColumn, ss[thisCopy.getColumnStatus(iColumn)],
995	primalSol[iColumn], primalDj[iColumn]);
996	#endif
997	// position at bound information
998	int jColumn = numberRows_;
999	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1000	double objValue = optimizationDirection_ * objective[iColumn];
1001	Status status = dualProblem->getRowStatus(iColumn);
1002	double otherValue = COIN_DBL_MAX;
1003	if (columnUpper_[iColumn] < 1.0e20 &&
1004	columnLower_[iColumn] > -1.0e20) {
1005	if (fabs(columnLower_[iColumn]) < fabs(columnUpper_[iColumn])) {
1006	otherValue = columnUpper_[iColumn] + dualDj[jColumn];
1007	} else {
1008	otherValue = columnLower_[iColumn] + dualDj[jColumn];
1009	}
1010	jColumn++;
1011	}
1012	if (status == basic) {
1013	// column is at bound
1014	if (otherValue == COIN_DBL_MAX) {
1015	reducedCost_[iColumn] = objValue - dualActs[iColumn];
1016	if (columnUpper_[iColumn] > 1.0e20) {
1017	if (columnLower_[iColumn] > -1.0e20) {
1018	if (columnUpper_[iColumn] > columnLower_[iColumn])
1019	setColumnStatus(iColumn, atLowerBound);
1020	else
1021	setColumnStatus(iColumn, isFixed);
1022	columnActivity_[iColumn] = columnLower_[iColumn];
1023	} else {
1024	// free
1025	setColumnStatus(iColumn, isFree);
1026	columnActivity_[iColumn] = 0.0;
1027	}
1028	} else {
1029	setColumnStatus(iColumn, atUpperBound);
1030	columnActivity_[iColumn] = columnUpper_[iColumn];
1031	}
1032	} else {
1033	reducedCost_[iColumn] = objValue - dualActs[iColumn];
1034	//printf("other dual sol %g\n",otherValue);
1035	if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1036	if (columnUpper_[iColumn] > columnLower_[iColumn])
1037	setColumnStatus(iColumn, atLowerBound);
1038	else
1039	setColumnStatus(iColumn, isFixed);
1040	columnActivity_[iColumn] = columnLower_[iColumn];
1041	} else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1042	if (columnUpper_[iColumn] > columnLower_[iColumn])
1043	setColumnStatus(iColumn, atUpperBound);
1044	else
1045	setColumnStatus(iColumn, isFixed);
1046	columnActivity_[iColumn] = columnUpper_[iColumn];
1047	} else {
1048	abort();
1049	}
1050	}
1051	} else {
1052	if (otherValue == COIN_DBL_MAX) {
1053	// column basic
1054	setColumnStatus(iColumn, basic);
1055	numberBasic++;
1056	if (columnLower_[iColumn] > -1.0e20) {
1057	columnActivity_[iColumn] = -dualDual[iColumn] + columnLower_[iColumn];
1058	} else if (columnUpper_[iColumn] < 1.0e20) {
1059	columnActivity_[iColumn] = -dualDual[iColumn] + columnUpper_[iColumn];
1060	} else {
1061	columnActivity_[iColumn] = -dualDual[iColumn];
1062	}
1063	reducedCost_[iColumn] = 0.0;
1064	} else {
1065	// may be at other bound
1066	//printf("xx %d %g jcol %d\n",iColumn,otherValue,jColumn-1);
1067	if (dualProblem->getColumnStatus(jColumn - 1) != basic) {
1068	// column basic
1069	setColumnStatus(iColumn, basic);
1070	numberBasic++;
1071	//printf("Col %d otherV %g dualDual %g\n",iColumn,
1072	// otherValue,dualDual[iColumn]);
1073	columnActivity_[iColumn] = -dualDual[iColumn];
1074	columnActivity_[iColumn] = otherValue;
1075	reducedCost_[iColumn] = 0.0;
1076	} else {
1077	reducedCost_[iColumn] = objValue - dualActs[iColumn];
1078	if (fabs(otherValue - columnLower_[iColumn]) < 1.0e-5) {
1079	if (columnUpper_[iColumn] > columnLower_[iColumn])
1080	setColumnStatus(iColumn, atLowerBound);
1081	else
1082	setColumnStatus(iColumn, isFixed);
1083	columnActivity_[iColumn] = columnLower_[iColumn];
1084	} else if (fabs(otherValue - columnUpper_[iColumn]) < 1.0e-5) {
1085	if (columnUpper_[iColumn] > columnLower_[iColumn])
1086	setColumnStatus(iColumn, atUpperBound);
1087	else
1088	setColumnStatus(iColumn, isFixed);
1089	columnActivity_[iColumn] = columnUpper_[iColumn];
1090	} else {
1091	abort();
1092	}
1093	}
1094	}
1095	}
1096	}
1097	// now rows
1098	int kExtraRow = jColumn;
1099	int numberRanges = 0;
1100	for (iRow = 0; iRow < numberRows_; iRow++) {
1101	Status status = dualProblem->getColumnStatus(iRow);
1102	if (status == basic) {
1103	// row is at bound
1104	dual_[iRow] = dualSol[iRow];;
1105	} else {
1106	// row basic
1107	setRowStatus(iRow, basic);
1108	numberBasic++;
1109	dual_[iRow] = 0.0;
1110	}
1111	if (rowLower_[iRow] < -1.0e20) {
1112	if (status == basic) {
1113	rowActivity_[iRow] = rowUpper_[iRow];
1114	setRowStatus(iRow, atUpperBound);
1115	} else {
1116	assert (dualDj[iRow] < 1.0e-5);
1117	rowActivity_[iRow] = rowUpper_[iRow] + dualDj[iRow];
1118	}
1119	} else if (rowUpper_[iRow] > 1.0e20) {
1120	if (status == basic) {
1121	rowActivity_[iRow] = rowLower_[iRow];
1122	setRowStatus(iRow, atLowerBound);
1123	} else {
1124	rowActivity_[iRow] = rowLower_[iRow] + dualDj[iRow];
1125	assert (dualDj[iRow] > -1.0e-5);
1126	}
1127	} else {
1128	if (rowUpper_[iRow] == rowLower_[iRow]) {
1129	rowActivity_[iRow] = rowLower_[iRow];
1130	if (status == basic) {
1131	setRowStatus(iRow, isFixed);
1132	}
1133	} else {
1134	// range
1135	numberRanges++;
1136	Status statusL = dualProblem->getColumnStatus(kExtraRow);
1137	//printf("range row %d (%d), extra %d (%d) - dualSol %g,%g dualDj %g,%g\n",
1138	// iRow,status,kExtraRow,statusL, dualSol[iRow],
1139	// dualSol[kExtraRow],dualDj[iRow],dualDj[kExtraRow]);
1140	if (status == basic) {
1141	assert (statusL != basic);
1142	rowActivity_[iRow] = rowUpper_[iRow];
1143	setRowStatus(iRow, atUpperBound);
1144	} else if (statusL == basic) {
1145	numberBasic--; // already counted
1146	rowActivity_[iRow] = rowLower_[iRow];
1147	setRowStatus(iRow, atLowerBound);
1148	dual_[iRow] = dualSol[kExtraRow];;
1149	} else {
1150	rowActivity_[iRow] = rowLower_[iRow] - dualDj[iRow];
1151	assert (dualDj[iRow] < 1.0e-5);
1152	// row basic
1153	//setRowStatus(iRow,basic);
1154	//numberBasic++;
1155	dual_[iRow] = 0.0;
1156	}
1157	kExtraRow++;
1158	}
1159	}
1160	}
1161	if (numberBasic != numberRows_) {
1162	printf("Bad basis - ranges - coding needed\n");
1163	assert (numberRanges);
1164	abort();
1165	}
1166	if (optimizationDirection_ < 0.0) {
1167	for (iRow = 0; iRow < numberRows_; iRow++) {
1168	dual_[iRow] = -dual_[iRow];
1169	}
1170	}
1171	// redo row activities
1172	memset(rowActivity_, 0, numberRows_ * sizeof(double));
1173	matrix_->times(1.0, columnActivity_, rowActivity_);
1174	// redo reduced costs
1175	memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
1176	matrix_->transposeTimes(-1.0, dual_, reducedCost_);
1177	checkSolutionInternal();
1178	if (sumDualInfeasibilities_ > 1.0e-5 \|\| sumPrimalInfeasibilities_ > 1.0e-5) {
1179	returnCode = 1;
1180	#ifdef CLP_INVESTIGATE
1181	printf("There are %d dual infeasibilities summing to %g ",
1182	numberDualInfeasibilities_, sumDualInfeasibilities_);
1183	printf("and %d primal infeasibilities summing to %g\n",
1184	numberPrimalInfeasibilities_, sumPrimalInfeasibilities_);
1185	#endif
1186	}
1187	// Below will go to ..DEBUG later
1188	#if 1 //ndef NDEBUG
1189	// Check if correct
1190	double * columnActivity = CoinCopyOfArray(columnActivity_, numberColumns_);
1191	double * rowActivity = CoinCopyOfArray(rowActivity_, numberRows_);
1192	double * reducedCost = CoinCopyOfArray(reducedCost_, numberColumns_);
1193	double * dual = CoinCopyOfArray(dual_, numberRows_);
1194	this->dual(); //primal();
1195	CoinRelFltEq eq(1.0e-5);
1196	for (iRow = 0; iRow < numberRows_; iRow++) {
1197	assert(eq(dual[iRow], dual_[iRow]));
1198	}
1199	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1200	assert(eq(columnActivity[iColumn], columnActivity_[iColumn]));
1201	}
1202	for (iRow = 0; iRow < numberRows_; iRow++) {
1203	assert(eq(rowActivity[iRow], rowActivity_[iRow]));
1204	}
1205	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1206	assert(eq(reducedCost[iColumn], reducedCost_[iColumn]));
1207	}
1208	delete [] columnActivity;
1209	delete [] rowActivity;
1210	delete [] reducedCost;
1211	delete [] dual;
1212	#endif
1213	return returnCode;
1214	}
1215	/* Does very cursory presolve.
1216	rhs is numberRows, whichRows is 3numberRows and whichColumns is 2numberColumns
1217	*/
1218	ClpSimplex *
1219	ClpSimplexOther::crunch(double * rhs, int * whichRow, int * whichColumn,
1220	int & nBound, bool moreBounds, bool tightenBounds)
1221	{
1222	//#define CHECK_STATUS
1223	#ifdef CHECK_STATUS
1224	{
1225	int n = 0;
1226	int i;
1227	for (i = 0; i < numberColumns_; i++)
1228	if (getColumnStatus(i) == ClpSimplex::basic)
1229	n++;
1230	for (i = 0; i < numberRows_; i++)
1231	if (getRowStatus(i) == ClpSimplex::basic)
1232	n++;
1233	assert (n == numberRows_);
1234	}
1235	#endif
1236
1237	const double * element = matrix_->getElements();
1238	const int * row = matrix_->getIndices();
1239	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1240	const int * columnLength = matrix_->getVectorLengths();
1241
1242	CoinZeroN(rhs, numberRows_);
1243	int iColumn;
1244	int iRow;
1245	CoinZeroN(whichRow, numberRows_);
1246	int * backColumn = whichColumn + numberColumns_;
1247	int numberRows2 = 0;
1248	int numberColumns2 = 0;
1249	double offset = 0.0;
1250	const double * objective = this->objective();
1251	double * solution = columnActivity_;
1252	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1253	double lower = columnLower_[iColumn];
1254	double upper = columnUpper_[iColumn];
1255	if (upper > lower \|\| getColumnStatus(iColumn) == ClpSimplex::basic) {
1256	backColumn[iColumn] = numberColumns2;
1257	whichColumn[numberColumns2++] = iColumn;
1258	for (CoinBigIndex j = columnStart[iColumn];
1259	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1260	int iRow = row[j];
1261	int n = whichRow[iRow];
1262	if (n == 0 && element[j])
1263	whichRow[iRow] = -iColumn - 1;
1264	else if (n < 0)
1265	whichRow[iRow] = 2;
1266	}
1267	} else {
1268	// fixed
1269	backColumn[iColumn] = -1;
1270	solution[iColumn] = upper;
1271	if (upper) {
1272	offset += objective[iColumn] * upper;
1273	for (CoinBigIndex j = columnStart[iColumn];
1274	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1275	int iRow = row[j];
1276	double value = element[j];
1277	rhs[iRow] += upper * value;
1278	}
1279	}
1280	}
1281	}
1282	int returnCode = 0;
1283	double tolerance = primalTolerance();
1284	nBound = 2 * numberRows_;
1285	for (iRow = 0; iRow < numberRows_; iRow++) {
1286	int n = whichRow[iRow];
1287	if (n > 0) {
1288	whichRow[numberRows2++] = iRow;
1289	} else if (n < 0) {
1290	//whichRow[numberRows2++]=iRow;
1291	//continue;
1292	// Can only do in certain circumstances as we don't know current value
1293	if (rowLower_[iRow] == rowUpper_[iRow] \|\| getRowStatus(iRow) == ClpSimplex::basic) {
1294	// save row and column for bound
1295	whichRow[--nBound] = iRow;
1296	whichRow[nBound+numberRows_] = -n - 1;
1297	} else if (moreBounds) {
1298	// save row and column for bound
1299	whichRow[--nBound] = iRow;
1300	whichRow[nBound+numberRows_] = -n - 1;
1301	} else {
1302	whichRow[numberRows2++] = iRow;
1303	}
1304	} else {
1305	// empty
1306	double rhsValue = rhs[iRow];
1307	if (rhsValue < rowLower_[iRow] - tolerance \|\| rhsValue > rowUpper_[iRow] + tolerance) {
1308	returnCode = 1; // infeasible
1309	}
1310	}
1311	}
1312	ClpSimplex * small = NULL;
1313	if (!returnCode) {
1314	//printf("CRUNCH from (%d,%d) to (%d,%d)\n",
1315	// numberRows_,numberColumns_,numberRows2,numberColumns2);
1316	small = new ClpSimplex(this, numberRows2, whichRow,
1317	numberColumns2, whichColumn, true, false);
1318	#if 0
1319	ClpPackedMatrix * rowCopy = dynamic_cast<ClpPackedMatrix *>(rowCopy_);
1320	if (rowCopy) {
1321	assert(!small->rowCopy());
1322	small->setNewRowCopy(new ClpPackedMatrix(*rowCopy, numberRows2, whichRow,
1323	numberColumns2, whichColumn));
1324	}
1325	#endif
1326	// Set some stuff
1327	small->setDualBound(dualBound_);
1328	small->setInfeasibilityCost(infeasibilityCost_);
1329	small->setSpecialOptions(specialOptions_);
1330	small->setPerturbation(perturbation_);
1331	small->defaultFactorizationFrequency();
1332	small->setAlphaAccuracy(alphaAccuracy_);
1333	// If no rows left then no tightening!
1334	if (!numberRows2 \|\| !numberColumns2)
1335	tightenBounds = false;
1336
1337	int numberElements = getNumElements();
1338	int numberElements2 = small->getNumElements();
1339	small->setObjectiveOffset(objectiveOffset() - offset);
1340	handler_->message(CLP_CRUNCH_STATS, messages_)
1341	<< numberRows2 << -(numberRows_ - numberRows2)
1342	<< numberColumns2 << -(numberColumns_ - numberColumns2)
1343	<< numberElements2 << -(numberElements - numberElements2)
1344	<< CoinMessageEol;
1345	// And set objective value to match
1346	small->setObjectiveValue(this->objectiveValue());
1347	double * rowLower2 = small->rowLower();
1348	double * rowUpper2 = small->rowUpper();
1349	int jRow;
1350	for (jRow = 0; jRow < numberRows2; jRow++) {
1351	iRow = whichRow[jRow];
1352	if (rowLower2[jRow] > -1.0e20)
1353	rowLower2[jRow] -= rhs[iRow];
1354	if (rowUpper2[jRow] < 1.0e20)
1355	rowUpper2[jRow] -= rhs[iRow];
1356	}
1357	// and bounds
1358	double * columnLower2 = small->columnLower();
1359	double * columnUpper2 = small->columnUpper();
1360	const char * integerInformation = integerType_;
1361	for (jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1362	iRow = whichRow[jRow];
1363	iColumn = whichRow[jRow+numberRows_];
1364	double lowerRow = rowLower_[iRow];
1365	if (lowerRow > -1.0e20)
1366	lowerRow -= rhs[iRow];
1367	double upperRow = rowUpper_[iRow];
1368	if (upperRow < 1.0e20)
1369	upperRow -= rhs[iRow];
1370	int jColumn = backColumn[iColumn];
1371	double lower = columnLower2[jColumn];
1372	double upper = columnUpper2[jColumn];
1373	double value = 0.0;
1374	for (CoinBigIndex j = columnStart[iColumn];
1375	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1376	if (iRow == row[j]) {
1377	value = element[j];
1378	break;
1379	}
1380	}
1381	assert (value);
1382	// convert rowLower and Upper to implied bounds on column
1383	double newLower = -COIN_DBL_MAX;
1384	double newUpper = COIN_DBL_MAX;
1385	if (value > 0.0) {
1386	if (lowerRow > -1.0e20)
1387	newLower = lowerRow / value;
1388	if (upperRow < 1.0e20)
1389	newUpper = upperRow / value;
1390	} else {
1391	if (upperRow < 1.0e20)
1392	newLower = upperRow / value;
1393	if (lowerRow > -1.0e20)
1394	newUpper = lowerRow / value;
1395	}
1396	if (integerInformation && integerInformation[iColumn]) {
1397	if (newLower - floor(newLower) < 10.0 * tolerance)
1398	newLower = floor(newLower);
1399	else
1400	newLower = ceil(newLower);
1401	if (ceil(newUpper) - newUpper < 10.0 * tolerance)
1402	newUpper = ceil(newUpper);
1403	else
1404	newUpper = floor(newUpper);
1405	}
1406	newLower = CoinMax(lower, newLower);
1407	newUpper = CoinMin(upper, newUpper);
1408	if (newLower > newUpper + tolerance) {
1409	//printf("XXYY inf on bound\n");
1410	returnCode = 1;
1411	}
1412	columnLower2[jColumn] = newLower;
1413	columnUpper2[jColumn] = CoinMax(newLower, newUpper);
1414	if (getRowStatus(iRow) != ClpSimplex::basic) {
1415	if (getColumnStatus(iColumn) == ClpSimplex::basic) {
1416	if (columnLower2[jColumn] == columnUpper2[jColumn]) {
1417	// can only get here if will be fixed
1418	small->setColumnStatus(jColumn, ClpSimplex::isFixed);
1419	} else {
1420	// solution is valid
1421	if (fabs(columnActivity_[iColumn] - columnLower2[jColumn]) <
1422	fabs(columnActivity_[iColumn] - columnUpper2[jColumn]))
1423	small->setColumnStatus(jColumn, ClpSimplex::atLowerBound);
1424	else
1425	small->setColumnStatus(jColumn, ClpSimplex::atUpperBound);
1426	}
1427	} else {
1428	//printf("what now neither basic\n");
1429	}
1430	}
1431	}
1432	if (returnCode) {
1433	delete small;
1434	small = NULL;
1435	} else if (tightenBounds && integerInformation) {
1436	// See if we can tighten any bounds
1437	// use rhs for upper and small duals for lower
1438	double * up = rhs;
1439	double * lo = small->dualRowSolution();
1440	const double * element = small->clpMatrix()->getElements();
1441	const int * row = small->clpMatrix()->getIndices();
1442	const CoinBigIndex * columnStart = small->clpMatrix()->getVectorStarts();
1443	//const int * columnLength = small->clpMatrix()->getVectorLengths();
1444	CoinZeroN(lo, numberRows2);
1445	CoinZeroN(up, numberRows2);
1446	for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1447	double upper = columnUpper2[iColumn];
1448	double lower = columnLower2[iColumn];
1449	//assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1450	for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1451	int iRow = row[j];
1452	double value = element[j];
1453	if (value > 0.0) {
1454	if (upper < 1.0e20)
1455	up[iRow] += upper * value;
1456	else
1457	up[iRow] = COIN_DBL_MAX;
1458	if (lower > -1.0e20)
1459	lo[iRow] += lower * value;
1460	else
1461	lo[iRow] = -COIN_DBL_MAX;
1462	} else {
1463	if (upper < 1.0e20)
1464	lo[iRow] += upper * value;
1465	else
1466	lo[iRow] = -COIN_DBL_MAX;
1467	if (lower > -1.0e20)
1468	up[iRow] += lower * value;
1469	else
1470	up[iRow] = COIN_DBL_MAX;
1471	}
1472	}
1473	}
1474	double * rowLower2 = small->rowLower();
1475	double * rowUpper2 = small->rowUpper();
1476	bool feasible = true;
1477	// make safer
1478	for (int iRow = 0; iRow < numberRows2; iRow++) {
1479	double lower = lo[iRow];
1480	if (lower > rowUpper2[iRow] + tolerance) {
1481	feasible = false;
1482	break;
1483	} else {
1484	lo[iRow] = CoinMin(lower - rowUpper2[iRow], 0.0) - tolerance;
1485	}
1486	double upper = up[iRow];
1487	if (upper < rowLower2[iRow] - tolerance) {
1488	feasible = false;
1489	break;
1490	} else {
1491	up[iRow] = CoinMax(upper - rowLower2[iRow], 0.0) + tolerance;
1492	}
1493	}
1494	if (!feasible) {
1495	delete small;
1496	small = NULL;
1497	} else {
1498	// and tighten
1499	for (int iColumn = 0; iColumn < numberColumns2; iColumn++) {
1500	if (integerInformation[whichColumn[iColumn]]) {
1501	double upper = columnUpper2[iColumn];
1502	double lower = columnLower2[iColumn];
1503	double newUpper = upper;
1504	double newLower = lower;
1505	double difference = upper - lower;
1506	if (lower > -1000.0 && upper < 1000.0) {
1507	for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1508	int iRow = row[j];
1509	double value = element[j];
1510	if (value > 0.0) {
1511	double upWithOut = up[iRow] - value * difference;
1512	if (upWithOut < 0.0) {
1513	newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1514	}
1515	double lowWithOut = lo[iRow] + value * difference;
1516	if (lowWithOut > 0.0) {
1517	newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1518	}
1519	} else {
1520	double upWithOut = up[iRow] + value * difference;
1521	if (upWithOut < 0.0) {
1522	newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1523	}
1524	double lowWithOut = lo[iRow] - value * difference;
1525	if (lowWithOut > 0.0) {
1526	newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1527	}
1528	}
1529	}
1530	if (newLower > lower \|\| newUpper < upper) {
1531	if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1532	newUpper = floor(newUpper);
1533	else
1534	newUpper = floor(newUpper + 0.5);
1535	if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1536	newLower = ceil(newLower);
1537	else
1538	newLower = ceil(newLower - 0.5);
1539	// change may be too small - check
1540	if (newLower > lower \|\| newUpper < upper) {
1541	if (newUpper >= newLower) {
1542	// Could also tighten in this
1543	//printf("%d bounds %g %g tightened to %g %g\n",
1544	// iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1545	// newLower,newUpper);
1546	#if 1
1547	columnUpper2[iColumn] = newUpper;
1548	columnLower2[iColumn] = newLower;
1549	columnUpper_[whichColumn[iColumn]] = newUpper;
1550	columnLower_[whichColumn[iColumn]] = newLower;
1551	#endif
1552	// and adjust bounds on rows
1553	newUpper -= upper;
1554	newLower -= lower;
1555	for (CoinBigIndex j = columnStart[iColumn]; j < columnStart[iColumn+1]; j++) {
1556	int iRow = row[j];
1557	double value = element[j];
1558	if (value > 0.0) {
1559	up[iRow] += newUpper * value;
1560	lo[iRow] += newLower * value;
1561	} else {
1562	lo[iRow] += newUpper * value;
1563	up[iRow] += newLower * value;
1564	}
1565	}
1566	} else {
1567	// infeasible
1568	//printf("%d bounds infeasible %g %g tightened to %g %g\n",
1569	// iColumn,columnLower2[iColumn],columnUpper2[iColumn],
1570	// newLower,newUpper);
1571	#if 1
1572	delete small;
1573	small = NULL;
1574	break;
1575	#endif
1576	}
1577	}
1578	}
1579	}
1580	}
1581	}
1582	}
1583	}
1584	}
1585	#if 0
1586	if (small) {
1587	static int which = 0;
1588	which++;
1589	char xxxx[20];
1590	sprintf(xxxx, "bad%d.mps", which);
1591	small->writeMps(xxxx, 0, 1);
1592	sprintf(xxxx, "largebad%d.mps", which);
1593	writeMps(xxxx, 0, 1);
1594	printf("bad%d %x old size %d %d new %d %d\n", which, small,
1595	numberRows_, numberColumns_, small->numberRows(), small->numberColumns());
1596	#if 0
1597	for (int i = 0; i < numberColumns_; i++)
1598	printf("Bound %d %g %g\n", i, columnLower_[i], columnUpper_[i]);
1599	for (int i = 0; i < numberRows_; i++)
1600	printf("Row bound %d %g %g\n", i, rowLower_[i], rowUpper_[i]);
1601	#endif
1602	}
1603	#endif
1604	#ifdef CHECK_STATUS
1605	{
1606	int n = 0;
1607	int i;
1608	for (i = 0; i < small->numberColumns(); i++)
1609	if (small->getColumnStatus(i) == ClpSimplex::basic)
1610	n++;
1611	for (i = 0; i < small->numberRows(); i++)
1612	if (small->getRowStatus(i) == ClpSimplex::basic)
1613	n++;
1614	assert (n == small->numberRows());
1615	}
1616	#endif
1617	return small;
1618	}
1619	/* After very cursory presolve.
1620	rhs is numberRows, whichRows is 3numberRows and whichColumns is 2numberColumns.
1621	*/
1622	void
1623	ClpSimplexOther::afterCrunch(const ClpSimplex & small,
1624	const int * whichRow,
1625	const int * whichColumn, int nBound)
1626	{
1627	#ifndef NDEBUG
1628	for (int i = 0; i < small.numberRows(); i++)
1629	assert (whichRow[i] >= 0 && whichRow[i] < numberRows_);
1630	for (int i = 0; i < small.numberColumns(); i++)
1631	assert (whichColumn[i] >= 0 && whichColumn[i] < numberColumns_);
1632	#endif
1633	getbackSolution(small, whichRow, whichColumn);
1634	// and deal with status for bounds
1635	const double * element = matrix_->getElements();
1636	const int * row = matrix_->getIndices();
1637	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1638	const int * columnLength = matrix_->getVectorLengths();
1639	double tolerance = primalTolerance();
1640	double djTolerance = dualTolerance();
1641	for (int jRow = nBound; jRow < 2 * numberRows_; jRow++) {
1642	int iRow = whichRow[jRow];
1643	int iColumn = whichRow[jRow+numberRows_];
1644	if (getColumnStatus(iColumn) != ClpSimplex::basic) {
1645	double lower = columnLower_[iColumn];
1646	double upper = columnUpper_[iColumn];
1647	double value = columnActivity_[iColumn];
1648	double djValue = reducedCost_[iColumn];
1649	dual_[iRow] = 0.0;
1650	if (upper > lower) {
1651	if (value < lower + tolerance && djValue > -djTolerance) {
1652	setColumnStatus(iColumn, ClpSimplex::atLowerBound);
1653	setRowStatus(iRow, ClpSimplex::basic);
1654	} else if (value > upper - tolerance && djValue < djTolerance) {
1655	setColumnStatus(iColumn, ClpSimplex::atUpperBound);
1656	setRowStatus(iRow, ClpSimplex::basic);
1657	} else {
1658	// has to be basic
1659	setColumnStatus(iColumn, ClpSimplex::basic);
1660	reducedCost_[iColumn] = 0.0;
1661	double value = 0.0;
1662	for (CoinBigIndex j = columnStart[iColumn];
1663	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1664	if (iRow == row[j]) {
1665	value = element[j];
1666	break;
1667	}
1668	}
1669	dual_[iRow] = djValue / value;
1670	if (rowUpper_[iRow] > rowLower_[iRow]) {
1671	if (fabs(rowActivity_[iRow] - rowLower_[iRow]) <
1672	fabs(rowActivity_[iRow] - rowUpper_[iRow]))
1673	setRowStatus(iRow, ClpSimplex::atLowerBound);
1674	else
1675	setRowStatus(iRow, ClpSimplex::atUpperBound);
1676	} else {
1677	setRowStatus(iRow, ClpSimplex::isFixed);
1678	}
1679	}
1680	} else {
1681	// row can always be basic
1682	setRowStatus(iRow, ClpSimplex::basic);
1683	}
1684	} else {
1685	// row can always be basic
1686	setRowStatus(iRow, ClpSimplex::basic);
1687	}
1688	}
1689	//#ifndef NDEBUG
1690	#if 0
1691	if (small.status() == 0) {
1692	int n = 0;
1693	int i;
1694	for (i = 0; i < numberColumns; i++)
1695	if (getColumnStatus(i) == ClpSimplex::basic)
1696	n++;
1697	for (i = 0; i < numberRows; i++)
1698	if (getRowStatus(i) == ClpSimplex::basic)
1699	n++;
1700	assert (n == numberRows);
1701	}
1702	#endif
1703	}
1704	/* Tightens integer bounds - returns number tightened or -1 if infeasible
1705	*/
1706	int
1707	ClpSimplexOther::tightenIntegerBounds(double * rhsSpace)
1708	{
1709	// See if we can tighten any bounds
1710	// use rhs for upper and small duals for lower
1711	double * up = rhsSpace;
1712	double * lo = dual_;
1713	const double * element = matrix_->getElements();
1714	const int * row = matrix_->getIndices();
1715	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
1716	const int * columnLength = matrix_->getVectorLengths();
1717	CoinZeroN(lo, numberRows_);
1718	CoinZeroN(up, numberRows_);
1719	for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1720	double upper = columnUpper_[iColumn];
1721	double lower = columnLower_[iColumn];
1722	//assert (columnLength[iColumn]==columnStart[iColumn+1]-columnStart[iColumn]);
1723	for (CoinBigIndex j = columnStart[iColumn];
1724	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1725	int iRow = row[j];
1726	double value = element[j];
1727	if (value > 0.0) {
1728	if (upper < 1.0e20)
1729	up[iRow] += upper * value;
1730	else
1731	up[iRow] = COIN_DBL_MAX;
1732	if (lower > -1.0e20)
1733	lo[iRow] += lower * value;
1734	else
1735	lo[iRow] = -COIN_DBL_MAX;
1736	} else {
1737	if (upper < 1.0e20)
1738	lo[iRow] += upper * value;
1739	else
1740	lo[iRow] = -COIN_DBL_MAX;
1741	if (lower > -1.0e20)
1742	up[iRow] += lower * value;
1743	else
1744	up[iRow] = COIN_DBL_MAX;
1745	}
1746	}
1747	}
1748	bool feasible = true;
1749	// make safer
1750	double tolerance = primalTolerance();
1751	for (int iRow = 0; iRow < numberRows_; iRow++) {
1752	double lower = lo[iRow];
1753	if (lower > rowUpper_[iRow] + tolerance) {
1754	feasible = false;
1755	break;
1756	} else {
1757	lo[iRow] = CoinMin(lower - rowUpper_[iRow], 0.0) - tolerance;
1758	}
1759	double upper = up[iRow];
1760	if (upper < rowLower_[iRow] - tolerance) {
1761	feasible = false;
1762	break;
1763	} else {
1764	up[iRow] = CoinMax(upper - rowLower_[iRow], 0.0) + tolerance;
1765	}
1766	}
1767	int numberTightened = 0;
1768	if (!feasible) {
1769	return -1;
1770	} else if (integerType_) {
1771	// and tighten
1772	for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
1773	if (integerType_[iColumn]) {
1774	double upper = columnUpper_[iColumn];
1775	double lower = columnLower_[iColumn];
1776	double newUpper = upper;
1777	double newLower = lower;
1778	double difference = upper - lower;
1779	if (lower > -1000.0 && upper < 1000.0) {
1780	for (CoinBigIndex j = columnStart[iColumn];
1781	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1782	int iRow = row[j];
1783	double value = element[j];
1784	if (value > 0.0) {
1785	double upWithOut = up[iRow] - value * difference;
1786	if (upWithOut < 0.0) {
1787	newLower = CoinMax(newLower, lower - (upWithOut + tolerance) / value);
1788	}
1789	double lowWithOut = lo[iRow] + value * difference;
1790	if (lowWithOut > 0.0) {
1791	newUpper = CoinMin(newUpper, upper - (lowWithOut - tolerance) / value);
1792	}
1793	} else {
1794	double upWithOut = up[iRow] + value * difference;
1795	if (upWithOut < 0.0) {
1796	newUpper = CoinMin(newUpper, upper - (upWithOut + tolerance) / value);
1797	}
1798	double lowWithOut = lo[iRow] - value * difference;
1799	if (lowWithOut > 0.0) {
1800	newLower = CoinMax(newLower, lower - (lowWithOut - tolerance) / value);
1801	}
1802	}
1803	}
1804	if (newLower > lower \|\| newUpper < upper) {
1805	if (fabs(newUpper - floor(newUpper + 0.5)) > 1.0e-6)
1806	newUpper = floor(newUpper);
1807	else
1808	newUpper = floor(newUpper + 0.5);
1809	if (fabs(newLower - ceil(newLower - 0.5)) > 1.0e-6)
1810	newLower = ceil(newLower);
1811	else
1812	newLower = ceil(newLower - 0.5);
1813	// change may be too small - check
1814	if (newLower > lower \|\| newUpper < upper) {
1815	if (newUpper >= newLower) {
1816	numberTightened++;
1817	//printf("%d bounds %g %g tightened to %g %g\n",
1818	// iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1819	// newLower,newUpper);
1820	columnUpper_[iColumn] = newUpper;
1821	columnLower_[iColumn] = newLower;
1822	// and adjust bounds on rows
1823	newUpper -= upper;
1824	newLower -= lower;
1825	for (CoinBigIndex j = columnStart[iColumn];
1826	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1827	int iRow = row[j];
1828	double value = element[j];
1829	if (value > 0.0) {
1830	up[iRow] += newUpper * value;
1831	lo[iRow] += newLower * value;
1832	} else {
1833	lo[iRow] += newUpper * value;
1834	up[iRow] += newLower * value;
1835	}
1836	}
1837	} else {
1838	// infeasible
1839	//printf("%d bounds infeasible %g %g tightened to %g %g\n",
1840	// iColumn,columnLower_[iColumn],columnUpper_[iColumn],
1841	// newLower,newUpper);
1842	return -1;
1843	}
1844	}
1845	}
1846	}
1847	}
1848	}
1849	}
1850	return numberTightened;
1851	}
1852	/* Parametrics
1853	This is an initial slow version.
1854	The code uses current bounds + theta * change (if change array not NULL)
1855	and similarly for objective.
1856	It starts at startingTheta and returns ending theta in endingTheta.
1857	If reportIncrement 0.0 it will report on any movement
1858	If reportIncrement >0.0 it will report at startingTheta+k*reportIncrement.
1859	If it can not reach input endingTheta return code will be 1 for infeasible,
1860	2 for unbounded, if error on ranges -1, otherwise 0.
1861	Normal report is just theta and objective but
1862	if event handler exists it may do more
1863	On exit endingTheta is maximum reached (can be used for next startingTheta)
1864	*/
1865	int
1866	ClpSimplexOther::parametrics(double startingTheta, double & endingTheta, double reportIncrement,
1867	const double * lowerChangeBound, const double * upperChangeBound,
1868	const double * lowerChangeRhs, const double * upperChangeRhs,
1869	const double * changeObjective)
1870	{
1871	bool needToDoSomething = true;
1872	bool canTryQuick = (reportIncrement) ? true : false;
1873	// Save copy of model
1874	ClpSimplex copyModel = *this;
1875	int savePerturbation = perturbation_;
1876	perturbation_ = 102; // switch off
1877	while (needToDoSomething) {
1878	needToDoSomething = false;
1879	algorithm_ = -1;
1880
1881	// save data
1882	ClpDataSave data = saveData();
1883	// Dantzig
1884	ClpDualRowPivot * savePivot = dualRowPivot_;
1885	dualRowPivot_ = new ClpDualRowDantzig();
1886	dualRowPivot_->setModel(this);
1887	int returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
1888	int iRow, iColumn;
1889	double * chgUpper = NULL;
1890	double * chgLower = NULL;
1891	double * chgObjective = NULL;
1892
1893
1894	if (!returnCode) {
1895	// Find theta when bounds will cross over and create arrays
1896	int numberTotal = numberRows_ + numberColumns_;
1897	chgLower = new double[numberTotal];
1898	memset(chgLower, 0, numberTotal * sizeof(double));
1899	chgUpper = new double[numberTotal];
1900	memset(chgUpper, 0, numberTotal * sizeof(double));
1901	chgObjective = new double[numberTotal];
1902	memset(chgObjective, 0, numberTotal * sizeof(double));
1903	assert (!rowScale_);
1904	double maxTheta = 1.0e50;
1905	if (lowerChangeRhs \|\| upperChangeRhs) {
1906	for (iRow = 0; iRow < numberRows_; iRow++) {
1907	double lower = rowLower_[iRow];
1908	double upper = rowUpper_[iRow];
1909	if (lower > upper) {
1910	maxTheta = -1.0;
1911	break;
1912	}
1913	double lowerChange = (lowerChangeRhs) ? lowerChangeRhs[iRow] : 0.0;
1914	double upperChange = (upperChangeRhs) ? upperChangeRhs[iRow] : 0.0;
1915	if (lower > -1.0e20 && upper < 1.0e20) {
1916	if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
1917	maxTheta = (upper - lower) / (lowerChange - upperChange);
1918	}
1919	}
1920	if (lower > -1.0e20) {
1921	lower_[numberColumns_+iRow] += startingTheta * lowerChange;
1922	chgLower[numberColumns_+iRow] = lowerChange;
1923	}
1924	if (upper < 1.0e20) {
1925	upper_[numberColumns_+iRow] += startingTheta * upperChange;
1926	chgUpper[numberColumns_+iRow] = upperChange;
1927	}
1928	}
1929	}
1930	if (maxTheta > 0.0) {
1931	if (lowerChangeBound \|\| upperChangeBound) {
1932	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1933	double lower = columnLower_[iColumn];
1934	double upper = columnUpper_[iColumn];
1935	if (lower > upper) {
1936	maxTheta = -1.0;
1937	break;
1938	}
1939	double lowerChange = (lowerChangeBound) ? lowerChangeBound[iColumn] : 0.0;
1940	double upperChange = (upperChangeBound) ? upperChangeBound[iColumn] : 0.0;
1941	if (lower > -1.0e20 && upper < 1.0e20) {
1942	if (lower + maxTheta * lowerChange > upper + maxTheta * upperChange) {
1943	maxTheta = (upper - lower) / (lowerChange - upperChange);
1944	}
1945	}
1946	if (lower > -1.0e20) {
1947	lower_[iColumn] += startingTheta * lowerChange;
1948	chgLower[iColumn] = lowerChange;
1949	}
1950	if (upper < 1.0e20) {
1951	upper_[iColumn] += startingTheta * upperChange;
1952	chgUpper[iColumn] = upperChange;
1953	}
1954	}
1955	}
1956	if (maxTheta == 1.0e50)
1957	maxTheta = COIN_DBL_MAX;
1958	}
1959	if (maxTheta < 0.0) {
1960	// bad ranges or initial
1961	returnCode = -1;
1962	}
1963	if (maxTheta < endingTheta) {
1964	char line[100];
1965	sprintf(line,"Crossover considerations reduce ending theta from %g to %g\n",
1966	endingTheta,maxTheta);
1967	handler_->message(CLP_GENERAL,messages_)
1968	<< line << CoinMessageEol;
1969	endingTheta = maxTheta;
1970	}
1971	if (endingTheta < startingTheta) {
1972	// bad initial
1973	returnCode = -2;
1974	}
1975	}
1976	double saveEndingTheta = endingTheta;
1977	if (!returnCode) {
1978	if (changeObjective) {
1979	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1980	chgObjective[iColumn] = changeObjective[iColumn];
1981	cost_[iColumn] += startingTheta * changeObjective[iColumn];
1982	}
1983	}
1984	double * saveDuals = NULL;
1985	reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
1986	assert (!problemStatus_);
1987	for (int i=0;i<numberRows_+numberColumns_;i++)
1988	setFakeBound(i, noFake);
1989	// Now do parametrics
1990	handler_->message(CLP_PARAMETRICS_STATS, messages_)
1991	<< startingTheta << objectiveValue() << CoinMessageEol;
1992	while (!returnCode) {
1993	//assert (reportIncrement);
1994	returnCode = parametricsLoop(startingTheta, endingTheta, reportIncrement,
1995	chgLower, chgUpper, chgObjective, data,
1996	canTryQuick);
1997	if (!returnCode) {
1998	//double change = endingTheta-startingTheta;
1999	startingTheta = endingTheta;
2000	endingTheta = saveEndingTheta;
2001	//for (int i=0;i<numberTotal;i++) {
2002	//lower_[i] += change*chgLower[i];
2003	//upper_[i] += change*chgUpper[i];
2004	//cost_[i] += change*chgObjective[i];
2005	//}
2006	handler_->message(CLP_PARAMETRICS_STATS, messages_)
2007	<< startingTheta << objectiveValue() << CoinMessageEol;
2008	if (startingTheta >= endingTheta)
2009	break;
2010	} else if (returnCode == -1) {
2011	// trouble - do external solve
2012	needToDoSomething = true;
2013	} else if (problemStatus_==1) {
2014	// can't move any further
2015	if (!canTryQuick) {
2016	handler_->message(CLP_PARAMETRICS_STATS, messages_)
2017	<< endingTheta << objectiveValue() << CoinMessageEol;
2018	problemStatus_=0;
2019	}
2020	} else {
2021	abort();
2022	}
2023	}
2024	}
2025	reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(0);
2026
2027	delete dualRowPivot_;
2028	dualRowPivot_ = savePivot;
2029	// Restore any saved stuff
2030	restoreData(data);
2031	if (needToDoSomething) {
2032	double saveStartingTheta = startingTheta; // known to be feasible
2033	int cleanedUp = 1;
2034	while (cleanedUp) {
2035	// tweak
2036	if (cleanedUp == 1) {
2037	if (!reportIncrement)
2038	startingTheta = CoinMin(startingTheta + 1.0e-5, saveEndingTheta);
2039	else
2040	startingTheta = CoinMin(startingTheta + reportIncrement, saveEndingTheta);
2041	} else {
2042	// restoring to go slowly
2043	startingTheta = saveStartingTheta;
2044	}
2045	// only works if not scaled
2046	int i;
2047	const double * obj1 = objective();
2048	double * obj2 = copyModel.objective();
2049	const double * lower1 = columnLower_;
2050	double * lower2 = copyModel.columnLower();
2051	const double * upper1 = columnUpper_;
2052	double * upper2 = copyModel.columnUpper();
2053	for (i = 0; i < numberColumns_; i++) {
2054	obj2[i] = obj1[i] + startingTheta * chgObjective[i];
2055	lower2[i] = lower1[i] + startingTheta * chgLower[i];
2056	upper2[i] = upper1[i] + startingTheta * chgUpper[i];
2057	}
2058	lower1 = rowLower_;
2059	lower2 = copyModel.rowLower();
2060	upper1 = rowUpper_;
2061	upper2 = copyModel.rowUpper();
2062	for (i = 0; i < numberRows_; i++) {
2063	lower2[i] = lower1[i] + startingTheta * chgLower[i+numberColumns_];
2064	upper2[i] = upper1[i] + startingTheta * chgUpper[i+numberColumns_];
2065	}
2066	copyModel.dual();
2067	if (copyModel.problemStatus()) {
2068	char line[100];
2069	sprintf(line,"Can not get to theta of %g\n", startingTheta);
2070	handler_->message(CLP_GENERAL,messages_)
2071	<< line << CoinMessageEol;
2072	canTryQuick = false; // do slowly to get exact amount
2073	// back to last known good
2074	if (cleanedUp == 1)
2075	cleanedUp = 2;
2076	else
2077	abort();
2078	} else {
2079	// and move stuff back
2080	int numberTotal = numberRows_ + numberColumns_;
2081	CoinMemcpyN(copyModel.statusArray(), numberTotal, status_);
2082	CoinMemcpyN(copyModel.primalColumnSolution(), numberColumns_, columnActivity_);
2083	CoinMemcpyN(copyModel.primalRowSolution(), numberRows_, rowActivity_);
2084	cleanedUp = 0;
2085	}
2086	}
2087	}
2088	delete [] chgLower;
2089	delete [] chgUpper;
2090	delete [] chgObjective;
2091	}
2092	perturbation_ = savePerturbation;
2093	char line[100];
2094	sprintf(line,"Ending theta %g\n", endingTheta);
2095	handler_->message(CLP_GENERAL,messages_)
2096	<< line << CoinMessageEol;
2097	return problemStatus_;
2098	}
2099	/* Version of parametrics which reads from file
2100	See CbcClpParam.cpp for details of format
2101	Returns -2 if unable to open file */
2102	int
2103	ClpSimplexOther::parametrics(const char * dataFile)
2104	{
2105	int returnCode=-2;
2106	FILE *fp = fopen(dataFile, "r");
2107	char line[200];
2108	if (!fp) {
2109	handler_->message(CLP_UNABLE_OPEN, messages_)
2110	<< dataFile << CoinMessageEol;
2111	return -2;
2112	}
2113
2114	if (!fgets(line, 200, fp)) {
2115	sprintf(line,"Empty parametrics file %s?",dataFile);
2116	handler_->message(CLP_GENERAL,messages_)
2117	<< line << CoinMessageEol;
2118	fclose(fp);
2119	return -2;
2120	}
2121	char * pos = line;
2122	char * put = line;
2123	while (pos >= ' ' && pos != '\n') {
2124	if (pos != ' ' && pos != '\t') {
2125	put = static_cast<char>(tolower(pos));
2126	put++;
2127	}
2128	pos++;
2129	}
2130	*put = '\0';
2131	pos = line;
2132	double startTheta=0.0;
2133	double endTheta=0.0;
2134	double intervalTheta=COIN_DBL_MAX;
2135	int detail=0;
2136	bool good = true;
2137	while (good) {
2138	good=false;
2139	// check ROWS
2140	char * comma = strchr(pos, ',');
2141	if (!comma)
2142	break;
2143	*comma = '\0';
2144	if (strcmp(pos,"rows"))
2145	break;
2146	*comma = ',';
2147	pos = comma+1;
2148	// check lower theta
2149	comma = strchr(pos, ',');
2150	if (!comma)
2151	break;
2152	*comma = '\0';
2153	startTheta = atof(pos);
2154	*comma = ',';
2155	pos = comma+1;
2156	// check upper theta
2157	comma = strchr(pos, ',');
2158	good=true;
2159	if (comma)
2160	*comma = '\0';
2161	endTheta = atof(pos);
2162	if (comma) {
2163	*comma = ',';
2164	pos = comma+1;
2165	comma = strchr(pos, ',');
2166	if (comma)
2167	*comma = '\0';
2168	intervalTheta = atof(pos);
2169	if (comma) {
2170	*comma = ',';
2171	pos = comma+1;
2172	comma = strchr(pos, ',');
2173	if (comma)
2174	*comma = '\0';
2175	detail = atoi(pos);
2176	if (comma)
2177	*comma = ',';
2178	}
2179	}
2180	break;
2181	}
2182	if (good) {
2183	if (startTheta<0.0\|\|
2184	startTheta>endTheta\|\|
2185	intervalTheta<0.0)
2186	good=false;
2187	if (detail<0\|\|detail>1)
2188	good=false;
2189	}
2190	if (intervalTheta>=endTheta)
2191	intervalTheta=0.0;
2192	if (!good) {
2193	sprintf(line,"Odd first line %s on file %s?",line,dataFile);
2194	handler_->message(CLP_GENERAL,messages_)
2195	<< line << CoinMessageEol;
2196	fclose(fp);
2197	return -2;
2198	}
2199	if (!fgets(line, 200, fp)) {
2200	sprintf(line,"Not enough records on parametrics file %s?",dataFile);
2201	handler_->message(CLP_GENERAL,messages_)
2202	<< line << CoinMessageEol;
2203	fclose(fp);
2204	return -2;
2205	}
2206	double * lowerRowMove = NULL;
2207	double * upperRowMove = NULL;
2208	double * lowerColumnMove = NULL;
2209	double * upperColumnMove = NULL;
2210	double * objectiveMove = NULL;
2211	char saveLine[200];
2212	saveLine[0]='\0';
2213	std::string headingsRow[] = {"name", "number", "lower", "upper", "rhs"};
2214	int gotRow[] = { -1, -1, -1, -1, -1};
2215	int orderRow[5];
2216	assert(sizeof(gotRow) == sizeof(orderRow));
2217	int nAcross = 0;
2218	pos = line;
2219	put = line;
2220	while (pos >= ' ' && pos != '\n') {
2221	if (pos != ' ' && pos != '\t') {
2222	put = static_cast<char>(tolower(pos));
2223	put++;
2224	}
2225	pos++;
2226	}
2227	*put = '\0';
2228	pos = line;
2229	int i;
2230	good = true;
2231	if (strncmp(line,"column",6)) {
2232	while (pos) {
2233	char * comma = strchr(pos, ',');
2234	if (comma)
2235	*comma = '\0';
2236	for (i = 0; i < static_cast<int> (sizeof(gotRow) / sizeof(int)); i++) {
2237	if (headingsRow[i] == pos) {
2238	if (gotRow[i] < 0) {
2239	orderRow[nAcross] = i;
2240	gotRow[i] = nAcross++;
2241	} else {
2242	// duplicate
2243	good = false;
2244	}
2245	break;
2246	}
2247	}
2248	if (i == static_cast<int> (sizeof(gotRow) / sizeof(int)))
2249	good = false;
2250	if (comma) {
2251	*comma = ',';
2252	pos = comma + 1;
2253	} else {
2254	break;
2255	}
2256	}
2257	if (gotRow[0] < 0 && gotRow[1] < 0)
2258	good = false;
2259	if (gotRow[0] >= 0 && gotRow[1] >= 0)
2260	good = false;
2261	if (gotRow[0] >= 0 && !lengthNames())
2262	good = false;
2263	if (gotRow[4]<0) {
2264	if (gotRow[2] < 0 && gotRow[3] >= 0)
2265	good = false;
2266	else if (gotRow[3] < 0 && gotRow[2] >= 0)
2267	good = false;
2268	} else if (gotRow[2]>=0\|\|gotRow[3]>=0) {
2269	good = false;
2270	}
2271	if (good) {
2272	char ** rowNames = new char * [numberRows_];
2273	int iRow;
2274	for (iRow = 0; iRow < numberRows_; iRow++) {
2275	rowNames[iRow] =
2276	CoinStrdup(rowName(iRow).c_str());
2277	}
2278	lowerRowMove = new double [numberRows_];
2279	memset(lowerRowMove,0,numberRows_*sizeof(double));
2280	upperRowMove = new double [numberRows_];
2281	memset(upperRowMove,0,numberRows_*sizeof(double));
2282	int nLine = 0;
2283	int nBadLine = 0;
2284	int nBadName = 0;
2285	bool goodLine=false;
2286	while (fgets(line, 200, fp)) {
2287	goodLine=true;
2288	if (!strncmp(line, "ENDATA", 6)\|\|
2289	!strncmp(line, "COLUMN",6))
2290	break;
2291	goodLine=false;
2292	nLine++;
2293	iRow = -1;
2294	double upper = 0.0;
2295	double lower = 0.0;
2296	char * pos = line;
2297	char * put = line;
2298	while (pos >= ' ' && pos != '\n') {
2299	if (pos != ' ' && pos != '\t') {
2300	put = pos;
2301	put++;
2302	}
2303	pos++;
2304	}
2305	*put = '\0';
2306	pos = line;
2307	for (int i = 0; i < nAcross; i++) {
2308	char * comma = strchr(pos, ',');
2309	if (comma) {
2310	*comma = '\0';
2311	} else if (i < nAcross - 1) {
2312	nBadLine++;
2313	break;
2314	}
2315	switch (orderRow[i]) {
2316	// name
2317	case 0:
2318	// For large problems this could be slow
2319	for (iRow = 0; iRow < numberRows_; iRow++) {
2320	if (!strcmp(rowNames[iRow], pos))
2321	break;
2322	}
2323	if (iRow == numberRows_)
2324	iRow = -1;
2325	break;
2326	// number
2327	case 1:
2328	iRow = atoi(pos);
2329	if (iRow < 0 \|\| iRow >= numberRows_)
2330	iRow = -1;
2331	break;
2332	// lower
2333	case 2:
2334	upper = atof(pos);
2335	break;
2336	// upper
2337	case 3:
2338	lower = atof(pos);
2339	break;
2340	// rhs
2341	case 4:
2342	lower = atof(pos);
2343	upper = lower;
2344	break;
2345	}
2346	if (comma) {
2347	*comma = ',';
2348	pos = comma + 1;
2349	}
2350	}
2351	if (iRow >= 0) {
2352	if (rowLower_[iRow]>-1.0e20)
2353	lowerRowMove[iRow] = lower;
2354	else
2355	lowerRowMove[iRow]=0.0;
2356	if (rowUpper_[iRow]<1.0e20)
2357	upperRowMove[iRow] = upper;
2358	else
2359	upperRowMove[iRow] = lower;
2360	} else {
2361	nBadName++;
2362	if(saveLine[0]=='\0')
2363	strcpy(saveLine,line);
2364	}
2365	}
2366	sprintf(line,"%d Row fields and %d records", nAcross, nLine);
2367	handler_->message(CLP_GENERAL,messages_)
2368	<< line << CoinMessageEol;
2369	if (nBadName) {
2370	sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2371	handler_->message(CLP_GENERAL,messages_)
2372	<< line << CoinMessageEol;
2373	returnCode=-1;
2374	good=false;
2375	}
2376	for (iRow = 0; iRow < numberRows_; iRow++) {
2377	free(rowNames[iRow]);
2378	}
2379	delete [] rowNames;
2380	} else {
2381	sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2382	handler_->message(CLP_GENERAL,messages_)
2383	<< line << CoinMessageEol;
2384	returnCode=-1;
2385	good=false;
2386	}
2387	}
2388	if (good&&(!strncmp(line, "COLUMN",6)\|\|!strncmp(line, "column",6))) {
2389	if (!fgets(line, 200, fp)) {
2390	sprintf(line,"Not enough records on parametrics file %s after COLUMNS?",dataFile);
2391	handler_->message(CLP_GENERAL,messages_)
2392	<< line << CoinMessageEol;
2393	fclose(fp);
2394	return -2;
2395	}
2396	std::string headingsColumn[] = {"name", "number", "lower", "upper", "objective"};
2397	saveLine[0]='\0';
2398	int gotColumn[] = { -1, -1, -1, -1, -1};
2399	int orderColumn[5];
2400	assert(sizeof(gotColumn) == sizeof(orderColumn));
2401	nAcross = 0;
2402	pos = line;
2403	put = line;
2404	while (pos >= ' ' && pos != '\n') {
2405	if (pos != ' ' && pos != '\t') {
2406	put = static_cast<char>(tolower(pos));
2407	put++;
2408	}
2409	pos++;
2410	}
2411	*put = '\0';
2412	pos = line;
2413	int i;
2414	if (strncmp(line,"endata",6)&&good) {
2415	while (pos) {
2416	char * comma = strchr(pos, ',');
2417	if (comma)
2418	*comma = '\0';
2419	for (i = 0; i < static_cast<int> (sizeof(gotColumn) / sizeof(int)); i++) {
2420	if (headingsColumn[i] == pos) {
2421	if (gotColumn[i] < 0) {
2422	orderColumn[nAcross] = i;
2423	gotColumn[i] = nAcross++;
2424	} else {
2425	// duplicate
2426	good = false;
2427	}
2428	break;
2429	}
2430	}
2431	if (i == static_cast<int> (sizeof(gotColumn) / sizeof(int)))
2432	good = false;
2433	if (comma) {
2434	*comma = ',';
2435	pos = comma + 1;
2436	} else {
2437	break;
2438	}
2439	}
2440	if (gotColumn[0] < 0 && gotColumn[1] < 0)
2441	good = false;
2442	if (gotColumn[0] >= 0 && gotColumn[1] >= 0)
2443	good = false;
2444	if (gotColumn[0] >= 0 && !lengthNames())
2445	good = false;
2446	if (good) {
2447	char ** columnNames = new char * [numberColumns_];
2448	int iColumn;
2449	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2450	columnNames[iColumn] =
2451	CoinStrdup(columnName(iColumn).c_str());
2452	}
2453	lowerColumnMove = reinterpret_cast<double > (malloc(numberColumns_ sizeof(double)));
2454	memset(lowerColumnMove,0,numberColumns_*sizeof(double));
2455	upperColumnMove = reinterpret_cast<double > (malloc(numberColumns_ sizeof(double)));
2456	memset(upperColumnMove,0,numberColumns_*sizeof(double));
2457	objectiveMove = reinterpret_cast<double > (malloc(numberColumns_ sizeof(double)));
2458	memset(objectiveMove,0,numberColumns_*sizeof(double));
2459	int nLine = 0;
2460	int nBadLine = 0;
2461	int nBadName = 0;
2462	bool goodLine=false;
2463	while (fgets(line, 200, fp)) {
2464	goodLine=true;
2465	if (!strncmp(line, "ENDATA", 6))
2466	break;
2467	goodLine=false;
2468	nLine++;
2469	iColumn = -1;
2470	double upper = 0.0;
2471	double lower = 0.0;
2472	double obj =0.0;
2473	char * pos = line;
2474	char * put = line;
2475	while (pos >= ' ' && pos != '\n') {
2476	if (pos != ' ' && pos != '\t') {
2477	put = pos;
2478	put++;
2479	}
2480	pos++;
2481	}
2482	*put = '\0';
2483	pos = line;
2484	for (int i = 0; i < nAcross; i++) {
2485	char * comma = strchr(pos, ',');
2486	if (comma) {
2487	*comma = '\0';
2488	} else if (i < nAcross - 1) {
2489	nBadLine++;
2490	break;
2491	}
2492	switch (orderColumn[i]) {
2493	// name
2494	case 0:
2495	// For large problems this could be slow
2496	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2497	if (!strcmp(columnNames[iColumn], pos))
2498	break;
2499	}
2500	if (iColumn == numberColumns_)
2501	iColumn = -1;
2502	break;
2503	// number
2504	case 1:
2505	iColumn = atoi(pos);
2506	if (iColumn < 0 \|\| iColumn >= numberColumns_)
2507	iColumn = -1;
2508	break;
2509	// lower
2510	case 2:
2511	upper = atof(pos);
2512	break;
2513	// upper
2514	case 3:
2515	lower = atof(pos);
2516	break;
2517	// objective
2518	case 4:
2519	obj = atof(pos);
2520	upper = lower;
2521	break;
2522	}
2523	if (comma) {
2524	*comma = ',';
2525	pos = comma + 1;
2526	}
2527	}
2528	if (iColumn >= 0) {
2529	if (columnLower_[iColumn]>-1.0e20)
2530	lowerColumnMove[iColumn] = lower;
2531	else
2532	lowerColumnMove[iColumn]=0.0;
2533	if (columnUpper_[iColumn]<1.0e20)
2534	upperColumnMove[iColumn] = upper;
2535	else
2536	upperColumnMove[iColumn] = lower;
2537	objectiveMove[iColumn] = obj;
2538	} else {
2539	nBadName++;
2540	if(saveLine[0]=='\0')
2541	strcpy(saveLine,line);
2542	}
2543	}
2544	sprintf(line,"%d Column fields and %d records", nAcross, nLine);
2545	handler_->message(CLP_GENERAL,messages_)
2546	<< line << CoinMessageEol;
2547	if (nBadName) {
2548	sprintf(line," ** %d records did not match on name/sequence, first bad %s", nBadName,saveLine);
2549	handler_->message(CLP_GENERAL,messages_)
2550	<< line << CoinMessageEol;
2551	returnCode=-1;
2552	good=false;
2553	}
2554	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2555	free(columnNames[iColumn]);
2556	}
2557	delete [] columnNames;
2558	} else {
2559	sprintf(line,"Duplicate or unknown keyword - or name/number fields wrong");
2560	handler_->message(CLP_GENERAL,messages_)
2561	<< line << CoinMessageEol;
2562	returnCode=-1;
2563	good=false;
2564	}
2565	}
2566	}
2567	returnCode=-1;
2568	if (good) {
2569	// clean arrays
2570	if (lowerRowMove) {
2571	bool empty=true;
2572	for (int i=0;i<numberRows_;i++) {
2573	if (lowerRowMove[i]) {
2574	empty=false;
2575	break;
2576	}
2577	}
2578	if (empty) {
2579	delete [] lowerRowMove;
2580	lowerRowMove=NULL;
2581	}
2582	}
2583	if (upperRowMove) {
2584	bool empty=true;
2585	for (int i=0;i<numberRows_;i++) {
2586	if (upperRowMove[i]) {
2587	empty=false;
2588	break;
2589	}
2590	}
2591	if (empty) {
2592	delete [] upperRowMove;
2593	upperRowMove=NULL;
2594	}
2595	}
2596	if (lowerColumnMove) {
2597	bool empty=true;
2598	for (int i=0;i<numberColumns_;i++) {
2599	if (lowerColumnMove[i]) {
2600	empty=false;
2601	break;
2602	}
2603	}
2604	if (empty) {
2605	delete [] lowerColumnMove;
2606	lowerColumnMove=NULL;
2607	}
2608	}
2609	if (upperColumnMove) {
2610	bool empty=true;
2611	for (int i=0;i<numberColumns_;i++) {
2612	if (upperColumnMove[i]) {
2613	empty=false;
2614	break;
2615	}
2616	}
2617	if (empty) {
2618	delete [] upperColumnMove;
2619	upperColumnMove=NULL;
2620	}
2621	}
2622	if (objectiveMove) {
2623	bool empty=true;
2624	for (int i=0;i<numberColumns_;i++) {
2625	if (objectiveMove[i]) {
2626	empty=false;
2627	break;
2628	}
2629	}
2630	if (empty) {
2631	delete [] objectiveMove;
2632	objectiveMove=NULL;
2633	}
2634	}
2635	int saveScaling = scalingFlag_;
2636	scalingFlag_ = 0;
2637	int saveLogLevel = handler_->logLevel();
2638	if (detail>0&&!intervalTheta)
2639	handler_->setLogLevel(3);
2640	else
2641	handler_->setLogLevel(1);
2642	returnCode = parametrics(startTheta,endTheta,intervalTheta,
2643	lowerColumnMove,upperColumnMove,
2644	lowerRowMove,upperRowMove,
2645	objectiveMove);
2646	scalingFlag_ = saveScaling;
2647	handler_->setLogLevel(saveLogLevel);
2648	}
2649	delete [] lowerRowMove;
2650	delete [] upperRowMove;
2651	delete [] lowerColumnMove;
2652	delete [] upperColumnMove;
2653	delete [] objectiveMove;
2654	fclose(fp);
2655	return returnCode;
2656	}
2657	int
2658	ClpSimplexOther::parametricsLoop(double startingTheta, double & endingTheta, double reportIncrement,
2659	const double * lowerChange, const double * upperChange,
2660	const double * changeObjective, ClpDataSave & data,
2661	bool canTryQuick)
2662	{
2663	// stuff is already at starting
2664	// For this crude version just try and go to end
2665	double change = 0.0;
2666	if (reportIncrement && canTryQuick) {
2667	endingTheta = CoinMin(endingTheta, startingTheta + reportIncrement);
2668	change = endingTheta - startingTheta;
2669	}
2670	int numberTotal = numberRows_ + numberColumns_;
2671	int i;
2672	for ( i = 0; i < numberTotal; i++) {
2673	lower_[i] += change * lowerChange[i];
2674	upper_[i] += change * upperChange[i];
2675	switch(getStatus(i)) {
2676
2677	case basic:
2678	case isFree:
2679	case superBasic:
2680	break;
2681	case isFixed:
2682	case atUpperBound:
2683	solution_[i] = upper_[i];
2684	break;
2685	case atLowerBound:
2686	solution_[i] = lower_[i];
2687	break;
2688	}
2689	cost_[i] += change * changeObjective[i];
2690	}
2691	problemStatus_ = -1;
2692
2693	// This says whether to restore things etc
2694	// startup will have factorized so can skip
2695	int factorType = 0;
2696	// Start check for cycles
2697	progress_.startCheck();
2698	// Say change made on first iteration
2699	changeMade_ = 1;
2700	/*
2701	Status of problem:
2702	0 - optimal
2703	1 - infeasible
2704	2 - unbounded
2705	-1 - iterating
2706	-2 - factorization wanted
2707	-3 - redo checking without factorization
2708	-4 - looks infeasible
2709	*/
2710	while (problemStatus_ < 0) {
2711	int iRow, iColumn;
2712	// clear
2713	for (iRow = 0; iRow < 4; iRow++) {
2714	rowArray_[iRow]->clear();
2715	}
2716
2717	for (iColumn = 0; iColumn < 2; iColumn++) {
2718	columnArray_[iColumn]->clear();
2719	}
2720
2721	// give matrix (and model costs and bounds a chance to be
2722	// refreshed (normally null)
2723	matrix_->refresh(this);
2724	// may factorize, checks if problem finished
2725	statusOfProblemInParametrics(factorType, data);
2726	// Say good factorization
2727	factorType = 1;
2728	if (data.sparseThreshold_) {
2729	// use default at present
2730	factorization_->sparseThreshold(0);
2731	factorization_->goSparse();
2732	}
2733
2734	// exit if victory declared
2735	if (problemStatus_ >= 0 &&
2736	(canTryQuick \|\| startingTheta>=endingTheta-1.0e-7) )
2737	break;
2738
2739	// test for maximum iterations
2740	if (hitMaximumIterations()) {
2741	problemStatus_ = 3;
2742	break;
2743	}
2744	// Check event
2745	{
2746	int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
2747	if (status >= 0) {
2748	problemStatus_ = 5;
2749	secondaryStatus_ = ClpEventHandler::endOfFactorization;
2750	break;
2751	}
2752	}
2753	// Do iterations
2754	problemStatus_=-1;
2755	if (canTryQuick) {
2756	double * saveDuals = NULL;
2757	reinterpret_cast<ClpSimplexDual *> (this)->whileIterating(saveDuals, 0);
2758	} else {
2759	whileIterating(startingTheta, endingTheta, reportIncrement,
2760	lowerChange, upperChange,
2761	changeObjective);
2762	startingTheta = endingTheta;
2763	}
2764	}
2765	if (!problemStatus_) {
2766	theta_ = change + startingTheta;
2767	eventHandler_->event(ClpEventHandler::theta);
2768	return 0;
2769	} else if (problemStatus_ == 10) {
2770	return -1;
2771	} else {
2772	return problemStatus_;
2773	}
2774	}
2775	/* Parametrics
2776	The code uses current bounds + theta * change (if change array not NULL)
2777	It starts at startingTheta and returns ending theta in endingTheta.
2778	If it can not reach input endingTheta return code will be 1 for infeasible,
2779	2 for unbounded, if error on ranges -1, otherwise 0.
2780	Event handler may do more
2781	On exit endingTheta is maximum reached (can be used for next startingTheta)
2782	*/
2783	int
2784	ClpSimplexOther::parametrics(double startingTheta, double & endingTheta,
2785	const double * lowerChangeBound, const double * upperChangeBound,
2786	const double * lowerChangeRhs, const double * upperChangeRhs)
2787	{
2788	int savePerturbation = perturbation_;
2789	perturbation_ = 102; // switch off
2790	algorithm_ = -1;
2791	// extra region
2792	int maximumPivots = factorization_->maximumPivots();
2793	int numberDense = factorization_->numberDense();
2794	int length = numberRows_ + numberDense + maximumPivots;
2795	assert (!rowArray_[4]);
2796	rowArray_[4]=new CoinIndexedVector(length);
2797	assert (!rowArray_[5]);
2798	rowArray_[5]=new CoinIndexedVector(length);
2799
2800	// save data
2801	ClpDataSave data = saveData();
2802	int numberTotal = numberRows_ + numberColumns_;
2803	int ratio = (2*sizeof(int))/sizeof(double);
2804	assert (ratio==1\|\|ratio==2);
2805	// allow for unscaled - even if not needed
2806	int lengthArrays = 4numberTotal+(numberTotal+2)ratio;
2807	/*
2808	Save information and modify
2809	*/
2810	double * saveLower = new double [lengthArrays];
2811	double * saveUpper = new double [lengthArrays];
2812	double * lowerCopy = saveLower+2*numberTotal;
2813	double * upperCopy = saveUpper+2*numberTotal;
2814	double * lowerChange = saveLower+numberTotal;
2815	double * upperChange = saveUpper+numberTotal;
2816	int * lowerList = (reinterpret_cast<int >(saveLower+4numberTotal))+2;
2817	int * upperList = (reinterpret_cast<int >(saveUpper+4numberTotal))+2;
2818	// To mark as odd
2819	char * markDone = reinterpret_cast<char *>(lowerList+numberTotal);
2820	memset(markDone,0,numberTotal);
2821	// Find theta when bounds will cross over and create arrays
2822	memset(lowerChange, 0, numberTotal * sizeof(double));
2823	memset(upperChange, 0, numberTotal * sizeof(double));
2824	if (lowerChangeBound)
2825	memcpy(lowerChange,lowerChangeBound,numberColumns_*sizeof(double));
2826	if (upperChangeBound)
2827	memcpy(upperChange,upperChangeBound,numberColumns_*sizeof(double));
2828	if (lowerChangeRhs)
2829	memcpy(lowerChange+numberColumns_,
2830	lowerChangeRhs,numberRows_*sizeof(double));
2831	if (upperChangeRhs)
2832	memcpy(upperChange+numberColumns_,
2833	upperChangeRhs,numberRows_*sizeof(double));
2834	int nLowerChange=0;
2835	int nUpperChange=0;
2836	for (int i=0;i<numberColumns_;i++) {
2837	if (lowerChange[i]) {
2838	lowerList[nLowerChange++]=i;
2839	}
2840	if (upperChange[i]) {
2841	upperList[nUpperChange++]=i;
2842	}
2843	}
2844	lowerList[-2]=nLowerChange;
2845	upperList[-2]=nUpperChange;
2846	for (int i=numberColumns_;i<numberTotal;i++) {
2847	if (lowerChange[i]) {
2848	lowerList[nLowerChange++]=i;
2849	}
2850	if (upperChange[i]) {
2851	upperList[nUpperChange++]=i;
2852	}
2853	}
2854	lowerList[-1]=nLowerChange;
2855	upperList[-1]=nUpperChange;
2856	memcpy(lowerCopy,columnLower_,numberColumns_*sizeof(double));
2857	memcpy(upperCopy,columnUpper_,numberColumns_*sizeof(double));
2858	memcpy(lowerCopy+numberColumns_,
2859	rowLower_,numberRows_*sizeof(double));
2860	memcpy(upperCopy+numberColumns_,
2861	rowUpper_,numberRows_*sizeof(double));
2862	{
2863	// extra for unscaled
2864	double * unscaledCopy;
2865	unscaledCopy = lowerCopy + numberTotal;
2866	memcpy(unscaledCopy,columnLower_,numberColumns_*sizeof(double));
2867	memcpy(unscaledCopy+numberColumns_,
2868	rowLower_,numberRows_*sizeof(double));
2869	unscaledCopy = upperCopy + numberTotal;
2870	memcpy(unscaledCopy,columnUpper_,numberColumns_*sizeof(double));
2871	memcpy(unscaledCopy+numberColumns_,
2872	rowUpper_,numberRows_*sizeof(double));
2873	}
2874	double maxTheta = 1.0e50;
2875	for (int iRow = 0; iRow < numberRows_; iRow++) {
2876	double lower = rowLower_[iRow];
2877	double upper = rowUpper_[iRow];
2878	if (lower<-1.0e30)
2879	lowerChange[numberColumns_+iRow]=0.0;
2880	double chgLower = lowerChange[numberColumns_+iRow];
2881	if (upper>1.0e30)
2882	upperChange[numberColumns_+iRow]=0.0;
2883	double chgUpper = upperChange[numberColumns_+iRow];
2884	if (lower > -1.0e30 && upper < 1.0e30) {
2885	if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2886	maxTheta = (upper - lower) / (chgLower - chgUpper);
2887	}
2888	}
2889	lower+=startingTheta*chgLower;
2890	upper+=startingTheta*chgUpper;
2891	if (lower > upper) {
2892	maxTheta = -1.0;
2893	break;
2894	}
2895	rowLower_[iRow]=lower;
2896	rowUpper_[iRow]=upper;
2897	}
2898	for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2899	double lower = columnLower_[iColumn];
2900	double upper = columnUpper_[iColumn];
2901	if (lower<-1.0e30)
2902	lowerChange[iColumn]=0.0;
2903	double chgLower = lowerChange[iColumn];
2904	if (upper>1.0e30)
2905	upperChange[iColumn]=0.0;
2906	double chgUpper = upperChange[iColumn];
2907	if (lower > -1.0e30 && upper < 1.0e30) {
2908	if (lower + maxTheta * chgLower > upper + maxTheta * chgUpper) {
2909	maxTheta = (upper - lower) / (chgLower - chgUpper);
2910	}
2911	}
2912	lower+=startingTheta*chgLower;
2913	upper+=startingTheta*chgUpper;
2914	if (lower > upper) {
2915	maxTheta = -1.0;
2916	break;
2917	}
2918	columnLower_[iColumn]=lower;
2919	columnUpper_[iColumn]=upper;
2920	}
2921	if (maxTheta == 1.0e50)
2922	maxTheta = COIN_DBL_MAX;
2923	int returnCode=0;
2924	if (maxTheta < 0.0) {
2925	// bad ranges or initial
2926	returnCode = -1;
2927	}
2928	if (maxTheta < endingTheta) {
2929	char line[100];
2930	sprintf(line,"Crossover considerations reduce ending theta from %g to %g\n",
2931	endingTheta,maxTheta);
2932	handler_->message(CLP_GENERAL,messages_)
2933	<< line << CoinMessageEol;
2934	endingTheta = maxTheta;
2935	}
2936	if (endingTheta < startingTheta) {
2937	// bad initial
2938	returnCode = -2;
2939	}
2940	bool swapped=false;
2941	// Dantzig
2942	#define ALL_DANTZIG
2943	#ifdef ALL_DANTZIG
2944	ClpDualRowPivot * savePivot = dualRowPivot_;
2945	dualRowPivot_ = new ClpDualRowDantzig();
2946	dualRowPivot_->setModel(this);
2947	#else
2948	ClpDualRowPivot * savePivot = NULL;
2949	#endif
2950	if (!returnCode) {
2951	returnCode = reinterpret_cast<ClpSimplexDual *> (this)->startupSolve(0, NULL, 0);
2952	if (!returnCode) {
2953	swapped=true;
2954	double * temp;
2955	memcpy(saveLower,lower_,numberTotal*sizeof(double));
2956	temp=saveLower;
2957	saveLower=lower_;
2958	lower_=temp;
2959	//columnLowerWork_ = lower_;
2960	//rowLowerWork_ = lower_ + numberColumns_;
2961	memcpy(saveUpper,upper_,numberTotal*sizeof(double));
2962	temp=saveUpper;
2963	saveUpper=upper_;
2964	upper_=temp;
2965	//columnUpperWork_ = upper_;
2966	//rowUpperWork_ = upper_ + numberColumns_;
2967	if (rowScale_) {
2968	// scale saved and change arrays
2969	double * lowerChange = lower_+numberTotal;
2970	double * upperChange = upper_+numberTotal;
2971	double * lowerSave = lowerChange+numberTotal;
2972	double * upperSave = upperChange+numberTotal;
2973	for (int i=0;i<numberColumns_;i++) {
2974	double multiplier = inverseColumnScale_[i];
2975	if (lowerSave[i]>-1.0e20)
2976	lowerSave[i] *= multiplier;
2977	if (upperSave[i]<1.0e20)
2978	upperSave[i] *= multiplier;
2979	lowerChange[i] *= multiplier;
2980	upperChange[i] *= multiplier;
2981	}
2982	lowerChange += numberColumns_;
2983	upperChange += numberColumns_;
2984	lowerSave += numberColumns_;
2985	upperSave += numberColumns_;
2986	for (int i=0;i<numberRows_;i++) {
2987	double multiplier = rowScale_[i];
2988	if (lowerSave[i]>-1.0e20)
2989	lowerSave[i] *= multiplier;
2990	if (upperSave[i]<1.0e20)
2991	upperSave[i] *= multiplier;
2992	lowerChange[i] *= multiplier;
2993	upperChange[i] *= multiplier;
2994	}
2995	}
2996	//double saveEndingTheta = endingTheta;
2997	double * saveDuals = NULL;
2998	reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
2999	if (numberPrimalInfeasibilities_&&sumPrimalInfeasibilities_<1.0e-4) {
3000	// probably can get rid of this if we adjust every change in theta
3001	//printf("INFEAS_A %d %g\n",numberPrimalInfeasibilities_,
3002	// sumPrimalInfeasibilities_);
3003	int pass=100;
3004	while(sumPrimalInfeasibilities_) {
3005	pass--;
3006	if (!pass)
3007	break;
3008	problemStatus_=-1;
3009	for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3010	double value=solution_[iSequence];
3011	// remember scaling
3012	if (value<lower_[iSequence]-1.0e-9) {
3013	lower_[iSequence]=value;
3014	lowerCopy[iSequence]=value;
3015	} else if (value>upper_[iSequence]+1.0e-9) {
3016	upper_[iSequence]=value;
3017	upperCopy[iSequence]=value;
3018	}
3019	}
3020	reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3021	}
3022	}
3023	assert (!problemStatus_);
3024	if (nLowerChange\|\|nUpperChange) {
3025	#ifndef ALL_DANTZIG
3026	// Dantzig
3027	savePivot = dualRowPivot_;
3028	dualRowPivot_ = new ClpDualRowDantzig();
3029	dualRowPivot_->setModel(this);
3030	#endif
3031	for (int i=0;i<numberRows_+numberColumns_;i++)
3032	setFakeBound(i, noFake);
3033	// Now do parametrics
3034	handler_->message(CLP_PARAMETRICS_STATS, messages_)
3035	<< startingTheta << objectiveValue() << CoinMessageEol;
3036	bool canSkipFactorization=true;
3037	while (!returnCode) {
3038	returnCode = parametricsLoop(startingTheta, endingTheta,
3039	data,canSkipFactorization);
3040	canSkipFactorization=false;
3041	if (!returnCode) {
3042	//startingTheta = endingTheta;
3043	//endingTheta = saveEndingTheta;
3044	handler_->message(CLP_PARAMETRICS_STATS, messages_)
3045	<< startingTheta << objectiveValue() << CoinMessageEol;
3046	if (startingTheta >= endingTheta-primalTolerance_
3047	\|\|problemStatus_==2)
3048	break;
3049	} else if (returnCode == -1) {
3050	// trouble - do external solve
3051	abort(); //needToDoSomething = true;
3052	} else if (problemStatus_==1) {
3053	// can't move any further
3054	handler_->message(CLP_PARAMETRICS_STATS, messages_)
3055	<< endingTheta << objectiveValue() << CoinMessageEol;
3056	problemStatus_=0;
3057	}
3058	}
3059	}
3060	//reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3061	}
3062	if (problemStatus_==2) {
3063	delete [] ray_;
3064	ray_ = new double [numberColumns_];
3065	}
3066	reinterpret_cast<ClpSimplexDual *> (this)->finishSolve(1);
3067	if (swapped&&lower_) {
3068	double * temp=saveLower;
3069	saveLower=lower_;
3070	lower_=temp;
3071	temp=saveUpper;
3072	saveUpper=upper_;
3073	upper_=temp;
3074	}
3075	}
3076	if (!scalingFlag_) {
3077	memcpy(columnLower_,lowerCopy,numberColumns_*sizeof(double));
3078	memcpy(columnUpper_,upperCopy,numberColumns_*sizeof(double));
3079	memcpy(rowLower_,lowerCopy+numberColumns_,
3080	numberRows_*sizeof(double));
3081	memcpy(rowUpper_,upperCopy+numberColumns_,
3082	numberRows_*sizeof(double));
3083	} else {
3084	// extra for unscaled
3085	double * unscaledCopy;
3086	unscaledCopy = lowerCopy + numberTotal;
3087	memcpy(columnLower_,unscaledCopy,numberColumns_*sizeof(double));
3088	memcpy(rowLower_,unscaledCopy+numberColumns_,
3089	numberRows_*sizeof(double));
3090	unscaledCopy = upperCopy + numberTotal;
3091	memcpy(columnUpper_,unscaledCopy,numberColumns_*sizeof(double));
3092	memcpy(rowUpper_,unscaledCopy+numberColumns_,
3093	numberRows_*sizeof(double));
3094	}
3095	delete [] saveLower;
3096	delete [] saveUpper;
3097	#ifdef ALL_DANTZIG
3098	if (savePivot) {
3099	#endif
3100	delete dualRowPivot_;
3101	dualRowPivot_ = savePivot;
3102	#ifdef ALL_DANTZIG
3103	}
3104	#endif
3105	// Restore any saved stuff
3106	restoreData(data);
3107	perturbation_ = savePerturbation;
3108	delete rowArray_[4];
3109	rowArray_[4]=NULL;
3110	delete rowArray_[5];
3111	rowArray_[5]=NULL;
3112	char line[100];
3113	sprintf(line,"Ending theta %g\n", endingTheta);
3114	handler_->message(CLP_GENERAL,messages_)
3115	<< line << CoinMessageEol;
3116	return problemStatus_;
3117	}
3118	int
3119	ClpSimplexOther::parametricsLoop(double & startingTheta, double & endingTheta,
3120	ClpDataSave & data,bool canSkipFactorization)
3121	{
3122	int numberTotal = numberRows_+numberColumns_;
3123	const double * lowerChange = lower_+numberTotal;
3124	const double * upperChange = upper_+numberTotal;
3125	// stuff is already at starting
3126	int * lowerList = (reinterpret_cast<int >(lower_+4numberTotal))+2;
3127	int * upperList = (reinterpret_cast<int >(upper_+4numberTotal))+2;
3128	problemStatus_ = -1;
3129	//double saveEndingTheta=endingTheta;
3130
3131	// This says whether to restore things etc
3132	// startup will have factorized so can skip
3133	int factorType = 0;
3134	// Start check for cycles
3135	progress_.startCheck();
3136	// Say change made on first iteration
3137	changeMade_ = 1;
3138	/*
3139	Status of problem:
3140	0 - optimal
3141	1 - infeasible
3142	2 - unbounded
3143	-1 - iterating
3144	-2 - factorization wanted
3145	-3 - redo checking without factorization
3146	-4 - looks infeasible
3147	*/
3148	while (problemStatus_ < 0) {
3149	int iRow, iColumn;
3150	// clear
3151	for (iRow = 0; iRow < 6; iRow++) {
3152	rowArray_[iRow]->clear();
3153	}
3154
3155	for (iColumn = 0; iColumn < 2; iColumn++) {
3156	columnArray_[iColumn]->clear();
3157	}
3158
3159	// give matrix (and model costs and bounds a chance to be
3160	// refreshed (normally null)
3161	matrix_->refresh(this);
3162	// may factorize, checks if problem finished
3163	if (!canSkipFactorization)
3164	statusOfProblemInParametrics(factorType, data);
3165	canSkipFactorization=false;
3166	if (numberPrimalInfeasibilities_) {
3167	// probably can get rid of this if we adjust every change in theta
3168	//printf("INFEAS %d %g\n",numberPrimalInfeasibilities_,
3169	// sumPrimalInfeasibilities_);
3170	const double * lowerChange = lower_+numberTotal;
3171	const double * upperChange = upper_+numberTotal;
3172	const double * startLower = lowerChange+numberTotal;
3173	const double * startUpper = upperChange+numberTotal;
3174	//startingTheta -= 1.0e-7;
3175	int nLowerChange = lowerList[-1];
3176	for (int i = 0; i < nLowerChange; i++) {
3177	int iSequence = lowerList[i];
3178	lower_[iSequence] = startLower[iSequence] + startingTheta * lowerChange[iSequence];
3179	}
3180	int nUpperChange = upperList[-1];
3181	for (int i = 0; i < nUpperChange; i++) {
3182	int iSequence = upperList[i];
3183	upper_[iSequence] = startUpper[iSequence] + startingTheta * upperChange[iSequence];
3184	}
3185	// adjust rhs in case dual uses
3186	memcpy(columnLower_,lower_,numberColumns_*sizeof(double));
3187	memcpy(rowLower_,lower_+numberColumns_,numberRows_*sizeof(double));
3188	memcpy(columnUpper_,upper_,numberColumns_*sizeof(double));
3189	memcpy(rowUpper_,upper_+numberColumns_,numberRows_*sizeof(double));
3190	if (rowScale_) {
3191	for (int i=0;i<numberColumns_;i++) {
3192	double multiplier = columnScale_[i];
3193	if (columnLower_[i]>-1.0e20)
3194	columnLower_[i] *= multiplier;
3195	if (columnUpper_[i]<1.0e20)
3196	columnUpper_[i] *= multiplier;
3197	}
3198	for (int i=0;i<numberRows_;i++) {
3199	double multiplier = inverseRowScale_[i];
3200	if (rowLower_[i]>-1.0e20)
3201	rowLower_[i] *= multiplier;
3202	if (rowUpper_[i]<1.0e20)
3203	rowUpper_[i] *= multiplier;
3204	}
3205	}
3206	double * saveDuals = NULL;
3207	problemStatus_=-1;
3208	reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(0, saveDuals, -1, data);
3209	int pass=100;
3210	double moved=0.0;
3211	while(sumPrimalInfeasibilities_) {
3212	pass--;
3213	if (!pass)
3214	break;
3215	problemStatus_=-1;
3216	for (int iSequence=numberColumns_;iSequence<numberTotal;iSequence++) {
3217	double value=solution_[iSequence];
3218	if (value<lower_[iSequence]-1.0e-9) {
3219	moved += lower_[iSequence]-value;
3220	lower_[iSequence]=value;
3221	} else if (value>upper_[iSequence]+1.0e-9) {
3222	moved = upper_[iSequence]-value;
3223	upper_[iSequence]=value;
3224	}
3225	}
3226	reinterpret_cast<ClpSimplexDual *> (this)->gutsOfDual(1, saveDuals, -1, data);
3227	}
3228	// adjust
3229	//printf("Should adjust - moved %g\n",moved);
3230	}
3231	// Say good factorization
3232	factorType = 1;
3233	if (data.sparseThreshold_) {
3234	// use default at present
3235	factorization_->sparseThreshold(0);
3236	factorization_->goSparse();
3237	}
3238
3239	// exit if victory declared
3240	if (problemStatus_ >= 0 && startingTheta>=endingTheta-1.0e-7 )
3241	break;
3242
3243	// test for maximum iterations
3244	if (hitMaximumIterations()) {
3245	problemStatus_ = 3;
3246	break;
3247	}
3248	#ifdef CLP_USER_DRIVEN
3249	// Check event
3250	{
3251	int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
3252	if (status >= 0) {
3253	problemStatus_ = 5;
3254	secondaryStatus_ = ClpEventHandler::endOfFactorization;
3255	break;
3256	}
3257	}
3258	#endif
3259	// Do iterations
3260	problemStatus_=-1;
3261	whileIterating(startingTheta, endingTheta, 0.0,
3262	lowerChange, upperChange,
3263	NULL);
3264	//startingTheta = endingTheta;
3265	//endingTheta = saveEndingTheta;
3266	}
3267	if (!problemStatus_/\|\|problemStatus_==2/) {
3268	theta_ = endingTheta;
3269	#ifdef CLP_USER_DRIVEN
3270	{
3271	double saveTheta=theta_;
3272	theta_ = endingTheta;
3273	int status=eventHandler_->event(ClpEventHandler::theta);
3274	if (status>=0&&status<10) {
3275	endingTheta=theta_;
3276	theta_=saveTheta;
3277	problemStatus_=-1;
3278	} else {
3279	if (status>=10) {
3280	problemStatus_=status-10;
3281	startingTheta=endingTheta;
3282	}
3283	theta_=saveTheta;
3284	}
3285	}
3286	#endif
3287	return 0;
3288	} else if (problemStatus_ == 10) {
3289	return -1;
3290	} else {
3291	return problemStatus_;
3292	}
3293	}
3294	/* Checks if finished. Updates status */
3295	void
3296	ClpSimplexOther::statusOfProblemInParametrics(int type, ClpDataSave & saveData)
3297	{
3298	if (type == 2) {
3299	// trouble - go to recovery
3300	problemStatus_ = 10;
3301	return;
3302	}
3303	if (problemStatus_ > -3 \|\| factorization_->pivots()) {
3304	// factorize
3305	// later on we will need to recover from singularities
3306	// also we could skip if first time
3307	if (type) {
3308	// is factorization okay?
3309	if (internalFactorize(1)) {
3310	// trouble - go to recovery
3311	problemStatus_ = 10;
3312	return;
3313	}
3314	}
3315	if (problemStatus_ != -4 \|\| factorization_->pivots() > 10)
3316	problemStatus_ = -3;
3317	}
3318	// at this stage status is -3 or -4 if looks infeasible
3319	// get primal and dual solutions
3320	gutsOfSolution(NULL, NULL);
3321	double realDualInfeasibilities = sumDualInfeasibilities_;
3322	// If bad accuracy treat as singular
3323	if ((largestPrimalError_ > 1.0e15 \|\| largestDualError_ > 1.0e15) && numberIterations_) {
3324	// trouble - go to recovery
3325	problemStatus_ = 10;
3326	return;
3327	} else if (largestPrimalError_ < 1.0e-7 && largestDualError_ < 1.0e-7) {
3328	// Can reduce tolerance
3329	double newTolerance = CoinMax(0.99 * factorization_->pivotTolerance(), saveData.pivotTolerance_);
3330	factorization_->pivotTolerance(newTolerance);
3331	}
3332	// Check if looping
3333	int loop;
3334	if (type != 2)
3335	loop = progress_.looping();
3336	else
3337	loop = -1;
3338	if (loop >= 0) {
3339	problemStatus_ = loop; //exit if in loop
3340	if (!problemStatus_) {
3341	// declaring victory
3342	numberPrimalInfeasibilities_ = 0;
3343	sumPrimalInfeasibilities_ = 0.0;
3344	} else {
3345	problemStatus_ = 10; // instead - try other algorithm
3346	}
3347	return;
3348	} else if (loop < -1) {
3349	// something may have changed
3350	gutsOfSolution(NULL, NULL);
3351	}
3352	progressFlag_ = 0; //reset progress flag
3353	if (handler_->detail(CLP_SIMPLEX_STATUS, messages_) < 100) {
3354	handler_->message(CLP_SIMPLEX_STATUS, messages_)
3355	<< numberIterations_ << objectiveValue();
3356	handler_->printing(sumPrimalInfeasibilities_ > 0.0)
3357	<< sumPrimalInfeasibilities_ << numberPrimalInfeasibilities_;
3358	handler_->printing(sumDualInfeasibilities_ > 0.0)
3359	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
3360	handler_->printing(numberDualInfeasibilitiesWithoutFree_
3361	< numberDualInfeasibilities_)
3362	<< numberDualInfeasibilitiesWithoutFree_;
3363	handler_->message() << CoinMessageEol;
3364	}
3365	/* If we are primal feasible and any dual infeasibilities are on
3366	free variables then it is better to go to primal */
3367	if (!numberPrimalInfeasibilities_ && !numberDualInfeasibilitiesWithoutFree_ &&
3368	numberDualInfeasibilities_) {
3369	problemStatus_ = 10;
3370	return;
3371	}
3372
3373	// check optimal
3374	// give code benefit of doubt
3375	if (sumOfRelaxedDualInfeasibilities_ == 0.0 &&
3376	sumOfRelaxedPrimalInfeasibilities_ == 0.0) {
3377	// say optimal (with these bounds etc)
3378	numberDualInfeasibilities_ = 0;
3379	sumDualInfeasibilities_ = 0.0;
3380	numberPrimalInfeasibilities_ = 0;
3381	sumPrimalInfeasibilities_ = 0.0;
3382	}
3383	if (dualFeasible() \|\| problemStatus_ == -4) {
3384	progress_.modifyObjective(objectiveValue_
3385	- sumDualInfeasibilities_ * dualBound_);
3386	}
3387	if (numberPrimalInfeasibilities_) {
3388	if (problemStatus_ == -4 \|\| problemStatus_ == -5) {
3389	problemStatus_ = 1; // infeasible
3390	}
3391	} else if (numberDualInfeasibilities_) {
3392	// clean up
3393	problemStatus_ = 10;
3394	} else {
3395	problemStatus_ = 0;
3396	}
3397	lastGoodIteration_ = numberIterations_;
3398	if (problemStatus_ < 0) {
3399	sumDualInfeasibilities_ = realDualInfeasibilities; // back to say be careful
3400	if (sumDualInfeasibilities_)
3401	numberDualInfeasibilities_ = 1;
3402	}
3403	// Allow matrices to be sorted etc
3404	int fake = -999; // signal sort
3405	matrix_->correctSequence(this, fake, fake);
3406	}
3407	//static double lastThetaX=0.0;
3408	/* This has the flow between re-factorizations
3409	Reasons to come out:
3410	-1 iterations etc
3411	-2 inaccuracy
3412	-3 slight inaccuracy (and done iterations)
3413	+0 looks optimal (might be unbounded - but we will investigate)
3414	+1 looks infeasible
3415	+3 max iterations
3416	+4 accuracy problems
3417	*/
3418	int
3419	ClpSimplexOther::whileIterating(double & startingTheta, double & endingTheta, double /reportIncrement/,
3420	const double * lowerChange, const double * upperChange,
3421	const double * /changeObjective/)
3422	{
3423	int numberTotal = numberColumns_ + numberRows_;
3424	const int * lowerList = (reinterpret_cast<int >(lower_+4numberTotal))+2;
3425	const int * upperList = (reinterpret_cast<int >(upper_+4numberTotal))+2;
3426	{
3427	int i;
3428	for (i = 0; i < 4; i++) {
3429	rowArray_[i]->clear();
3430	}
3431	for (i = 0; i < 2; i++) {
3432	columnArray_[i]->clear();
3433	}
3434	}
3435	// if can't trust much and long way from optimal then relax
3436	if (largestPrimalError_ > 10.0)
3437	factorization_->relaxAccuracyCheck(CoinMin(1.0e2, largestPrimalError_ / 10.0));
3438	else
3439	factorization_->relaxAccuracyCheck(1.0);
3440	// status stays at -1 while iterating, >=0 finished, -2 to invert
3441	// status -3 to go to top without an invert
3442	int returnCode = -1;
3443	double lastTheta = startingTheta;
3444	double useTheta = startingTheta;
3445	while (problemStatus_ == -1) {
3446	double increaseTheta = CoinMin(endingTheta - lastTheta, 1.0e50);
3447	// Get theta for bounds - we know can't crossover
3448	int pivotType = nextTheta(1, increaseTheta,
3449	lowerChange, upperChange, NULL);
3450	useTheta += theta_;
3451	double change = useTheta - lastTheta;
3452	if (change>1.0e-14) {
3453	int n;
3454	n=lowerList[-1];
3455	for (int i=0;i<n;i++) {
3456	int iSequence = lowerList[i];
3457	lower_[iSequence] += change * lowerChange[iSequence];
3458	if(getStatus(iSequence)==atLowerBound) {
3459	solution_[iSequence] = lower_[iSequence];
3460	}
3461	}
3462	n=upperList[-1];
3463	for (int i=0;i<n;i++) {
3464	int iSequence = upperList[i];
3465	upper_[iSequence] += change * upperChange[iSequence];
3466	if(getStatus(iSequence)==atUpperBound\|\|
3467	getStatus(iSequence)==isFixed) {
3468	solution_[iSequence] = upper_[iSequence];
3469	}
3470	}
3471	}
3472	sequenceIn_=-1;
3473	if (pivotType) {
3474	if (useTheta>lastTheta+1.0e-9) {
3475	handler_->message(CLP_PARAMETRICS_STATS, messages_)
3476	<< useTheta << objectiveValue() << CoinMessageEol;
3477	lastTheta = useTheta;
3478	}
3479	problemStatus_ = -2;
3480	if (!factorization_->pivots()&&pivotRow_<0)
3481	problemStatus_=2;
3482	#ifdef CLP_USER_DRIVEN
3483	{
3484	double saveTheta=theta_;
3485	theta_ = endingTheta;
3486	int status=eventHandler_->event(ClpEventHandler::theta);
3487	if (status>=0&&status<10) {
3488	endingTheta=theta_;
3489	theta_=saveTheta;
3490	problemStatus_=-1;
3491	continue;
3492	} else {
3493	if (status>=10)
3494	problemStatus_=status-10;
3495	if (status<0)
3496	startingTheta = useTheta;
3497	theta_=saveTheta;
3498	}
3499	}
3500	#else
3501	startingTheta = useTheta;
3502	#endif
3503	return 4;
3504	}
3505	// choose row to go out
3506	//reinterpret_cast<ClpSimplexDual *> ( this)->dualRow(-1);
3507	if (pivotRow_ >= 0) {
3508	// we found a pivot row
3509	if (handler_->detail(CLP_SIMPLEX_PIVOTROW, messages_) < 100) {
3510	handler_->message(CLP_SIMPLEX_PIVOTROW, messages_)
3511	<< pivotRow_
3512	<< CoinMessageEol;
3513	}
3514	// check accuracy of weights
3515	dualRowPivot_->checkAccuracy();
3516	// do ratio test for normal iteration
3517	double bestPossiblePivot = bestPivot();
3518	if (sequenceIn_ >= 0) {
3519	// normal iteration
3520	// update the incoming column
3521	double btranAlpha = -alpha_ * directionOut_; // for check
3522	unpackPacked(rowArray_[1]);
3523	// and update dual weights (can do in parallel - with extra array)
3524	rowArray_[2]->clear();
3525	alpha_ = dualRowPivot_->updateWeights(rowArray_[0],
3526	rowArray_[2],
3527	rowArray_[3],
3528	rowArray_[1]);
3529	// see if update stable
3530	#ifdef CLP_DEBUG
3531	if ((handler_->logLevel() & 32))
3532	printf("btran alpha %g, ftran alpha %g\n", btranAlpha, alpha_);
3533	#endif
3534	double checkValue = 1.0e-7;
3535	// if can't trust much and long way from optimal then relax
3536	if (largestPrimalError_ > 10.0)
3537	checkValue = CoinMin(1.0e-4, 1.0e-8 * largestPrimalError_);
3538	if (fabs(btranAlpha) < 1.0e-12 \|\| fabs(alpha_) < 1.0e-12 \|\|
3539	fabs(btranAlpha - alpha_) > checkValue*(1.0 + fabs(alpha_))) {
3540	handler_->message(CLP_DUAL_CHECK, messages_)
3541	<< btranAlpha
3542	<< alpha_
3543	<< CoinMessageEol;
3544	// clear arrays
3545	rowArray_[4]->clear();
3546	if (factorization_->pivots()) {
3547	dualRowPivot_->unrollWeights();
3548	problemStatus_ = -2; // factorize now
3549	rowArray_[0]->clear();
3550	rowArray_[1]->clear();
3551	columnArray_[0]->clear();
3552	returnCode = -2;
3553	break;
3554	} else {
3555	// take on more relaxed criterion
3556	double test;
3557	if (fabs(btranAlpha) < 1.0e-8 \|\| fabs(alpha_) < 1.0e-8)
3558	test = 1.0e-1 * fabs(alpha_);
3559	else
3560	test = 1.0e-4 * (1.0 + fabs(alpha_));
3561	if (fabs(btranAlpha) < 1.0e-12 \|\| fabs(alpha_) < 1.0e-12 \|\|
3562	fabs(btranAlpha - alpha_) > test) {
3563	dualRowPivot_->unrollWeights();
3564	// need to reject something
3565	char x = isColumn(sequenceOut_) ? 'C' : 'R';
3566	handler_->message(CLP_SIMPLEX_FLAG, messages_)
3567	<< x << sequenceWithin(sequenceOut_)
3568	<< CoinMessageEol;
3569	setFlagged(sequenceOut_);
3570	progress_.clearBadTimes();
3571	lastBadIteration_ = numberIterations_; // say be more cautious
3572	rowArray_[0]->clear();
3573	rowArray_[1]->clear();
3574	columnArray_[0]->clear();
3575	if (fabs(alpha_) < 1.0e-10 && fabs(btranAlpha) < 1.0e-8 && numberIterations_ > 100) {
3576	//printf("I think should declare infeasible\n");
3577	problemStatus_ = 1;
3578	returnCode = 1;
3579	break;
3580	}
3581	continue;
3582	}
3583	}
3584	}
3585	// update duals BEFORE replaceColumn so can do updateColumn
3586	double objectiveChange = 0.0;
3587	// do duals first as variables may flip bounds
3588	// rowArray_[0] and columnArray_[0] may have flips
3589	// so use rowArray_[3] for work array from here on
3590	int nswapped = 0;
3591	//rowArray_[0]->cleanAndPackSafe(1.0e-60);
3592	//columnArray_[0]->cleanAndPackSafe(1.0e-60);
3593	if ((specialOptions_&2097152)==0) {
3594	nswapped = reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0],
3595	rowArray_[2], theta_,
3596	objectiveChange, false);
3597	assert (!nswapped);
3598	} else {
3599	rowArray_[0]->clear();
3600	rowArray_[2]->clear();
3601	columnArray_[0]->clear();
3602	}
3603	// which will change basic solution
3604	if (nswapped) {
3605	factorization_->updateColumn(rowArray_[3], rowArray_[2]);
3606	dualRowPivot_->updatePrimalSolution(rowArray_[2],
3607	1.0, objectiveChange);
3608	// recompute dualOut_
3609	valueOut_ = solution_[sequenceOut_];
3610	if (directionOut_ < 0) {
3611	dualOut_ = valueOut_ - upperOut_;
3612	} else {
3613	dualOut_ = lowerOut_ - valueOut_;
3614	}
3615	}
3616	// amount primal will move
3617	double movement = -dualOut_ * directionOut_ / alpha_;
3618	// so objective should increase by fabs(dj)*movement
3619	// but we already have objective change - so check will be good
3620	if (objectiveChange + fabs(movement * dualIn_) < -1.0e-5) {
3621	#ifdef CLP_DEBUG
3622	if (handler_->logLevel() & 32)
3623	printf("movement %g, swap change %g, rest %g * %g\n",
3624	objectiveChange + fabs(movement * dualIn_),
3625	objectiveChange, movement, dualIn_);
3626	#endif
3627	assert (objectiveChange + fabs(movement * dualIn_) >= -1.0e-5);
3628	if(factorization_->pivots()) {
3629	// going backwards - factorize
3630	dualRowPivot_->unrollWeights();
3631	problemStatus_ = -2; // factorize now
3632	returnCode = -2;
3633	break;
3634	}
3635	}
3636	CoinAssert(fabs(dualOut_) < 1.0e50);
3637	// if stable replace in basis
3638	int updateStatus = factorization_->replaceColumn(this,
3639	rowArray_[2],
3640	rowArray_[1],
3641	pivotRow_,
3642	alpha_);
3643	// if no pivots, bad update but reasonable alpha - take and invert
3644	if (updateStatus == 2 &&
3645	!factorization_->pivots() && fabs(alpha_) > 1.0e-5)
3646	updateStatus = 4;
3647	if (updateStatus == 1 \|\| updateStatus == 4) {
3648	// slight error
3649	if (factorization_->pivots() > 5 \|\| updateStatus == 4) {
3650	problemStatus_ = -2; // factorize now
3651	returnCode = -3;
3652	}
3653	} else if (updateStatus == 2) {
3654	// major error
3655	dualRowPivot_->unrollWeights();
3656	// later we may need to unwind more e.g. fake bounds
3657	if (factorization_->pivots()) {
3658	problemStatus_ = -2; // factorize now
3659	returnCode = -2;
3660	break;
3661	} else {
3662	// need to reject something
3663	char x = isColumn(sequenceOut_) ? 'C' : 'R';
3664	handler_->message(CLP_SIMPLEX_FLAG, messages_)
3665	<< x << sequenceWithin(sequenceOut_)
3666	<< CoinMessageEol;
3667	setFlagged(sequenceOut_);
3668	progress_.clearBadTimes();
3669	lastBadIteration_ = numberIterations_; // say be more cautious
3670	rowArray_[0]->clear();
3671	rowArray_[1]->clear();
3672	columnArray_[0]->clear();
3673	// make sure dual feasible
3674	// look at all rows and columns
3675	double objectiveChange = 0.0;
3676	reinterpret_cast<ClpSimplexDual *> ( this)->updateDualsInDual(rowArray_[0], columnArray_[0], rowArray_[1],
3677	0.0, objectiveChange, true);
3678	continue;
3679	}
3680	} else if (updateStatus == 3) {
3681	// out of memory
3682	// increase space if not many iterations
3683	if (factorization_->pivots() <
3684	0.5 * factorization_->maximumPivots() &&
3685	factorization_->pivots() < 200)
3686	factorization_->areaFactor(
3687	factorization_->areaFactor() * 1.1);
3688	problemStatus_ = -2; // factorize now
3689	} else if (updateStatus == 5) {
3690	problemStatus_ = -2; // factorize now
3691	}
3692	// update change vector
3693	{
3694	double * work = rowArray_[1]->denseVector();
3695	int number = rowArray_[1]->getNumElements();
3696	int * which = rowArray_[1]->getIndices();
3697	assert (!rowArray_[4]->packedMode());
3698	assert (rowArray_[1]->packedMode());
3699	double pivotValue = rowArray_[4]->denseVector()[pivotRow_];
3700	double multiplier = -pivotValue/alpha_;
3701	if (multiplier) {
3702	for (int i = 0; i < number; i++) {
3703	int iRow = which[i];
3704	rowArray_[4]->quickAddNonZero(iRow,multiplier*work[i]);
3705	}
3706	}
3707	pivotValue = rowArray_[4]->denseVector()[pivotRow_];
3708	// we want pivot to be -multiplier
3709	rowArray_[4]->quickAdd(pivotRow_,-multiplier-pivotValue);
3710	}
3711	// update primal solution
3712	if (theta_ < 0.0) {
3713	#ifdef CLP_DEBUG
3714	if (handler_->logLevel() & 32)
3715	printf("negative theta %g\n", theta_);
3716	#endif
3717	theta_ = 0.0;
3718	}
3719	// do actual flips
3720	reinterpret_cast<ClpSimplexDual *> ( this)->flipBounds(rowArray_[0], columnArray_[0]);
3721	//rowArray_[1]->expand();
3722	dualRowPivot_->updatePrimalSolution(rowArray_[1],
3723	movement,
3724	objectiveChange);
3725	// modify dualout
3726	dualOut_ /= alpha_;
3727	dualOut_ *= -directionOut_;
3728	//setStatus(sequenceIn_,basic);
3729	dj_[sequenceIn_] = 0.0;
3730	//double oldValue = valueIn_;
3731	if (directionIn_ == -1) {
3732	// as if from upper bound
3733	valueIn_ = upperIn_ + dualOut_;
3734	} else {
3735	// as if from lower bound
3736	valueIn_ = lowerIn_ + dualOut_;
3737	}
3738	objectiveChange = 0.0;
3739	if ((specialOptions_&2097152)==0) {
3740	for (int i=0;i<numberTotal;i++)
3741	objectiveChange += solution_[i]*cost_[i];
3742	objectiveChange -= objectiveValue_;
3743	}
3744	// outgoing
3745	originalBound(sequenceOut_,useTheta,lowerChange,upperChange);
3746	lowerOut_=lower_[sequenceOut_];
3747	upperOut_=upper_[sequenceOut_];
3748	// set dj to zero unless values pass
3749	if (directionOut_ > 0) {
3750	valueOut_ = lowerOut_;
3751	dj_[sequenceOut_] = theta_;
3752	} else {
3753	valueOut_ = upperOut_;
3754	dj_[sequenceOut_] = -theta_;
3755	}
3756	solution_[sequenceOut_] = valueOut_;
3757	int whatNext = housekeeping(objectiveChange);
3758	{
3759	char in[200],out[200];
3760	int iSequence=sequenceIn_;
3761	if (iSequence<numberColumns_) {
3762	if (lengthNames_)
3763	strcpy(in,columnNames_[iSequence].c_str());
3764	else
3765	sprintf(in,"C%7.7d",iSequence);
3766	} else {
3767	iSequence -= numberColumns_;
3768	if (lengthNames_)
3769	strcpy(in,rowNames_[iSequence].c_str());
3770	else
3771	sprintf(in,"R%7.7d",iSequence);
3772	}
3773	iSequence=sequenceOut_;
3774	if (iSequence<numberColumns_) {
3775	if (lengthNames_)
3776	strcpy(out,columnNames_[iSequence].c_str());
3777	else
3778	sprintf(out,"C%7.7d",iSequence);
3779	} else {
3780	iSequence -= numberColumns_;
3781	if (lengthNames_)
3782	strcpy(out,rowNames_[iSequence].c_str());
3783	else
3784	sprintf(out,"R%7.7d",iSequence);
3785	}
3786	handler_->message(CLP_PARAMETRICS_STATS2, messages_)
3787	<< useTheta << objectiveValue()
3788	<< in << out << CoinMessageEol;
3789	}
3790	if (useTheta>lastTheta+1.0e-9) {
3791	handler_->message(CLP_PARAMETRICS_STATS, messages_)
3792	<< useTheta << objectiveValue() << CoinMessageEol;
3793	lastTheta = useTheta;
3794	}
3795	// and set bounds correctly
3796	originalBound(sequenceIn_,useTheta,lowerChange,upperChange);
3797	reinterpret_cast<ClpSimplexDual *> ( this)->changeBound(sequenceOut_);
3798	if (whatNext == 1) {
3799	problemStatus_ = -2; // refactorize
3800	} else if (whatNext == 2) {
3801	// maximum iterations or equivalent
3802	problemStatus_ = 3;
3803	returnCode = 3;
3804	break;
3805	}
3806	#ifdef CLP_USER_DRIVEN
3807	// Check event
3808	{
3809	int status = eventHandler_->event(ClpEventHandler::endOfIteration);
3810	if (status >= 0) {
3811	problemStatus_ = 5;
3812	secondaryStatus_ = ClpEventHandler::endOfIteration;
3813	returnCode = 4;
3814	break;
3815	}
3816	}
3817	#endif
3818	} else {
3819	// no incoming column is valid
3820	#ifdef CLP_USER_DRIVEN
3821	rowArray_[0]->clear();
3822	columnArray_[0]->clear();
3823	theta_ = useTheta;
3824	lastTheta = useTheta;
3825	int action = eventHandler_->event(ClpEventHandler::noTheta);
3826	if (action>=0) {
3827	endingTheta=theta_;
3828	theta_ = 0.0;
3829	//adjust [4] from handler - but
3830	rowArray_[4]->clear(); // temp
3831	if (action>=0&&action<10)
3832	problemStatus_=-1; // carry on
3833	else if (action==15)
3834	problemStatus_ =5; // say stopped
3835	returnCode = 1;
3836	if (action==0\|\|action>=10)
3837	break;
3838	else
3839	continue;
3840	} else {
3841	theta_ = 0.0;
3842	}
3843	#endif
3844	pivotRow_ = -1;
3845	#ifdef CLP_DEBUG
3846	if (handler_->logLevel() & 32)
3847	printf("** no column pivot\n");
3848	#endif
3849	if (factorization_->pivots() < 10) {
3850	// If we have just factorized and infeasibility reasonable say infeas
3851	if (((specialOptions_ & 4096) != 0 \|\| bestPossiblePivot < 1.0e-11) && dualBound_ > 1.0e8) {
3852	if (valueOut_ > upperOut_ + 1.0e-3 \|\| valueOut_ < lowerOut_ - 1.0e-3
3853	\|\| (specialOptions_ & 64) == 0) {
3854	// say infeasible
3855	problemStatus_ = 1;
3856	// unless primal feasible!!!!
3857	//printf("%d %g %d %g\n",numberPrimalInfeasibilities_,sumPrimalInfeasibilities_,
3858	// numberDualInfeasibilities_,sumDualInfeasibilities_);
3859	if (numberDualInfeasibilities_)
3860	problemStatus_ = 10;
3861	rowArray_[0]->clear();
3862	columnArray_[0]->clear();
3863	}
3864	}
3865	// If special option set - put off as long as possible
3866	if ((specialOptions_ & 64) == 0) {
3867	problemStatus_ = -4; //say looks infeasible
3868	} else {
3869	// flag
3870	char x = isColumn(sequenceOut_) ? 'C' : 'R';
3871	handler_->message(CLP_SIMPLEX_FLAG, messages_)
3872	<< x << sequenceWithin(sequenceOut_)
3873	<< CoinMessageEol;
3874	setFlagged(sequenceOut_);
3875	if (!factorization_->pivots()) {
3876	rowArray_[0]->clear();
3877	columnArray_[0]->clear();
3878	continue;
3879	}
3880	}
3881	}
3882	rowArray_[0]->clear();
3883	columnArray_[0]->clear();
3884	returnCode = 1;
3885	break;
3886	}
3887	} else {
3888	// no pivot row
3889	#ifdef CLP_USER_DRIVEN
3890	{
3891	double saveTheta=theta_;
3892	theta_ = endingTheta;
3893	int status=eventHandler_->event(ClpEventHandler::theta);
3894	if (status>=0&&status<10) {
3895	endingTheta=theta_;
3896	theta_=saveTheta;
3897	continue;
3898	} else {
3899	theta_=saveTheta;
3900	}
3901	}
3902	#endif
3903	#ifdef CLP_DEBUG
3904	if (handler_->logLevel() & 32)
3905	printf("** no row pivot\n");
3906	#endif
3907	int numberPivots = factorization_->pivots();
3908	bool specialCase;
3909	int useNumberFake;
3910	returnCode = 0;
3911	if (numberPivots < 20 &&
3912	(specialOptions_ & 2048) != 0 && !numberChanged_ && perturbation_ >= 100
3913	&& dualBound_ > 1.0e8) {
3914	specialCase = true;
3915	// as dual bound high - should be okay
3916	useNumberFake = 0;
3917	} else {
3918	specialCase = false;
3919	useNumberFake = numberFake_;
3920	}
3921	if (!numberPivots \|\| specialCase) {
3922	// may have crept through - so may be optimal
3923	// check any flagged variables
3924	int iRow;
3925	for (iRow = 0; iRow < numberRows_; iRow++) {
3926	int iPivot = pivotVariable_[iRow];
3927	if (flagged(iPivot))
3928	break;
3929	}
3930	if (iRow < numberRows_ && numberPivots) {
3931	// try factorization
3932	returnCode = -2;
3933	}
3934
3935	if (useNumberFake \|\| numberDualInfeasibilities_) {
3936	// may be dual infeasible
3937	problemStatus_ = -5;
3938	} else {
3939	if (iRow < numberRows_) {
3940	problemStatus_ = -5;
3941	} else {
3942	if (numberPivots) {
3943	// objective may be wrong
3944	objectiveValue_ = innerProduct(cost_,
3945	numberColumns_ + numberRows_,
3946	solution_);
3947	objectiveValue_ += objective_->nonlinearOffset();
3948	objectiveValue_ /= (objectiveScale_ * rhsScale_);
3949	if ((specialOptions_ & 16384) == 0) {
3950	// and dual_ may be wrong (i.e. for fixed or basic)
3951	CoinIndexedVector * arrayVector = rowArray_[1];
3952	arrayVector->clear();
3953	int iRow;
3954	double * array = arrayVector->denseVector();
3955	/* Use dual_ instead of array
3956	Even though dual_ is only numberRows_ long this is
3957	okay as gets permuted to longer rowArray_[2]
3958	*/
3959	arrayVector->setDenseVector(dual_);
3960	int * index = arrayVector->getIndices();
3961	int number = 0;
3962	for (iRow = 0; iRow < numberRows_; iRow++) {
3963	int iPivot = pivotVariable_[iRow];
3964	double value = cost_[iPivot];
3965	dual_[iRow] = value;
3966	if (value) {
3967	index[number++] = iRow;
3968	}
3969	}
3970	arrayVector->setNumElements(number);
3971	// Extended duals before "updateTranspose"
3972	matrix_->dualExpanded(this, arrayVector, NULL, 0);
3973	// Btran basic costs
3974	rowArray_[2]->clear();
3975	factorization_->updateColumnTranspose(rowArray_[2], arrayVector);
3976	// and return vector
3977	arrayVector->setDenseVector(array);
3978	}
3979	}
3980	problemStatus_ = 0;
3981	sumPrimalInfeasibilities_ = 0.0;
3982	if ((specialOptions_&(1024 + 16384)) != 0) {
3983	CoinIndexedVector * arrayVector = rowArray_[1];
3984	arrayVector->clear();
3985	double * rhs = arrayVector->denseVector();
3986	times(1.0, solution_, rhs);
3987	bool bad2 = false;
3988	int i;
3989	for ( i = 0; i < numberRows_; i++) {
3990	if (rhs[i] < rowLowerWork_[i] - primalTolerance_ \|\|
3991	rhs[i] > rowUpperWork_[i] + primalTolerance_) {
3992	bad2 = true;
3993	} else if (fabs(rhs[i] - rowActivityWork_[i]) > 1.0e-3) {
3994	}
3995	rhs[i] = 0.0;
3996	}
3997	for ( i = 0; i < numberColumns_; i++) {
3998	if (solution_[i] < columnLowerWork_[i] - primalTolerance_ \|\|
3999	solution_[i] > columnUpperWork_[i] + primalTolerance_) {
4000	bad2 = true;
4001	}
4002	}
4003	if (bad2) {
4004	problemStatus_ = -3;
4005	returnCode = -2;
4006	// Force to re-factorize early next time
4007	int numberPivots = factorization_->pivots();
4008	forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4009	}
4010	}
4011	}
4012	}
4013	} else {
4014	problemStatus_ = -3;
4015	returnCode = -2;
4016	// Force to re-factorize early next time
4017	int numberPivots = factorization_->pivots();
4018	forceFactorization_ = CoinMin(forceFactorization_, (numberPivots + 1) >> 1);
4019	}
4020	break;
4021	}
4022	}
4023	startingTheta = lastTheta+theta_;
4024	return returnCode;
4025	}
4026	// Finds best possible pivot
4027	double
4028	ClpSimplexOther::bestPivot(bool justColumns)
4029	{
4030	// Get good size for pivot
4031	// Allow first few iterations to take tiny
4032	double acceptablePivot = 1.0e-9;
4033	if (numberIterations_ > 100)
4034	acceptablePivot = 1.0e-8;
4035	if (factorization_->pivots() > 10 \|\|
4036	(factorization_->pivots() && sumDualInfeasibilities_))
4037	acceptablePivot = 1.0e-5; // if we have iterated be more strict
4038	else if (factorization_->pivots() > 5)
4039	acceptablePivot = 1.0e-6; // if we have iterated be slightly more strict
4040	else if (factorization_->pivots())
4041	acceptablePivot = 1.0e-8; // relax
4042	double bestPossiblePivot = 1.0;
4043	// get sign for finding row of tableau
4044	// normal iteration
4045	// create as packed
4046	double direction = directionOut_;
4047	rowArray_[0]->createPacked(1, &pivotRow_, &direction);
4048	factorization_->updateColumnTranspose(rowArray_[1], rowArray_[0]);
4049	// put row of tableau in rowArray[0] and columnArray[0]
4050	matrix_->transposeTimes(this, -1.0,
4051	rowArray_[0], rowArray_[3], columnArray_[0]);
4052	sequenceIn_=-1;
4053	if (justColumns)
4054	rowArray_[0]->clear();
4055	// do ratio test for normal iteration
4056	bestPossiblePivot =
4057	reinterpret_cast<ClpSimplexDual *>
4058	( this)->dualColumn(rowArray_[0],
4059	columnArray_[0], columnArray_[1],
4060	rowArray_[3], acceptablePivot, NULL);
4061	return bestPossiblePivot;
4062	}
4063	// Computes next theta and says if objective or bounds (0= bounds, 1 objective, -1 none)
4064	int
4065	ClpSimplexOther::nextTheta(int type, double maxTheta,
4066	const double * lowerChange, const double * upperChange,
4067	const double * /changeObjective/)
4068	{
4069	int numberTotal = numberColumns_ + numberRows_;
4070	const int * lowerList = (reinterpret_cast<int >(lower_+4numberTotal))+2;
4071	const int * upperList = (reinterpret_cast<int >(upper_+4numberTotal))+2;
4072	int iSequence;
4073	theta_ = maxTheta;
4074	bool toLower = false;
4075	assert (type==1);
4076	// may need to decide based on model?
4077	bool needFullUpdate = rowArray_[4]->getNumElements()==0;
4078	double * array = rowArray_[4]->denseVector();
4079	const int * row = matrix_->getIndices();
4080	const int * columnLength = matrix_->getVectorLengths();
4081	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4082	const double * elementByColumn = matrix_->getElements();
4083	if (!factorization_->pivots()\|\|needFullUpdate) {
4084	rowArray_[4]->clear();
4085	// get change
4086	if (!rowScale_) {
4087	int n;
4088	n=lowerList[-2];
4089	int i;
4090	for (i=0;i<n;i++) {
4091	int iSequence = lowerList[i];
4092	assert (iSequence<numberColumns_);
4093	if (getColumnStatus(iSequence)==atLowerBound) {
4094	double value=lowerChange[iSequence];
4095	for (CoinBigIndex j = columnStart[iSequence];
4096	j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4097	rowArray_[4]->quickAddNonZero(row[j], elementByColumn[j]*value);
4098	}
4099	}
4100	}
4101	n=lowerList[-1];
4102	const double * change = lowerChange+numberColumns_;
4103	for (;i<n;i++) {
4104	int iSequence = lowerList[i]-numberColumns_;
4105	assert (iSequence>=0);
4106	if (getRowStatus(iSequence)==atLowerBound) {
4107	double value=change[iSequence];
4108	rowArray_[4]->quickAddNonZero(iSequence, -value);
4109	}
4110	}
4111	n=upperList[-2];
4112	for (i=0;i<n;i++) {
4113	int iSequence = upperList[i];
4114	assert (iSequence<numberColumns_);
4115	if (getColumnStatus(iSequence)==atUpperBound) {
4116	double value=upperChange[iSequence];
4117	for (CoinBigIndex j = columnStart[iSequence];
4118	j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4119	rowArray_[4]->quickAddNonZero(row[j], elementByColumn[j]*value);
4120	}
4121	}
4122	}
4123	n=upperList[-1];
4124	change = upperChange+numberColumns_;
4125	for (;i<n;i++) {
4126	int iSequence = upperList[i]-numberColumns_;
4127	assert (iSequence>=0);
4128	if (getRowStatus(iSequence)==atUpperBound) {
4129	double value=change[iSequence];
4130	rowArray_[4]->quickAddNonZero(iSequence, -value);
4131	}
4132	}
4133	} else {
4134	int n;
4135	n=lowerList[-2];
4136	int i;
4137	for (i=0;i<n;i++) {
4138	int iSequence = lowerList[i];
4139	assert (iSequence<numberColumns_);
4140	if (getColumnStatus(iSequence)==atLowerBound) {
4141	double value=lowerChange[iSequence];
4142	// apply scaling
4143	double scale = columnScale_[iSequence];
4144	for (CoinBigIndex j = columnStart[iSequence];
4145	j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4146	int iRow = row[j];
4147	rowArray_[4]->quickAddNonZero(iRow, elementByColumn[j]scale rowScale_[iRow]*value);
4148	}
4149	}
4150	}
4151	n=lowerList[-1];
4152	const double * change = lowerChange+numberColumns_;
4153	for (;i<n;i++) {
4154	int iSequence = lowerList[i]-numberColumns_;
4155	assert (iSequence>=0);
4156	if (getRowStatus(iSequence)==atLowerBound) {
4157	double value=change[iSequence];
4158	rowArray_[4]->quickAddNonZero(iSequence, -value);
4159	}
4160	}
4161	n=upperList[-2];
4162	for (i=0;i<n;i++) {
4163	int iSequence = upperList[i];
4164	assert (iSequence<numberColumns_);
4165	if (getColumnStatus(iSequence)==atUpperBound) {
4166	double value=upperChange[iSequence];
4167	// apply scaling
4168	double scale = columnScale_[iSequence];
4169	for (CoinBigIndex j = columnStart[iSequence];
4170	j < columnStart[iSequence] + columnLength[iSequence]; j++) {
4171	int iRow = row[j];
4172	rowArray_[4]->quickAddNonZero(iRow, elementByColumn[j]scale rowScale_[iRow]*value);
4173	}
4174	}
4175	}
4176	n=upperList[-1];
4177	change = upperChange+numberColumns_;
4178	for (;i<n;i++) {
4179	int iSequence = upperList[i]-numberColumns_;
4180	assert (iSequence>=0);
4181	if (getRowStatus(iSequence)==atUpperBound) {
4182	double value=change[iSequence];
4183	rowArray_[4]->quickAddNonZero(iSequence, -value);
4184	}
4185	}
4186	}
4187	// ftran it
4188	factorization_->updateColumn(rowArray_[0], rowArray_[4]);
4189	} else {
4190	assert (sequenceIn_>=0);
4191	assert (sequenceOut_>=0);
4192	assert (sequenceIn_!=sequenceOut_);
4193	double change = (directionIn_>0) ? -lowerChange[sequenceIn_] : -upperChange[sequenceIn_];
4194	int needed=0;
4195	assert (!rowArray_[5]->getNumElements());
4196	if (change) {
4197	if (sequenceIn_<numberColumns_) {
4198	if (!rowScale_) {
4199	for (CoinBigIndex i = columnStart[sequenceIn_];
4200	i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4201	rowArray_[5]->quickAddNonZero(row[i], elementByColumn[i]*change);
4202	}
4203	} else {
4204	// apply scaling
4205	double scale = columnScale_[sequenceIn_];
4206	for (CoinBigIndex i = columnStart[sequenceIn_];
4207	i < columnStart[sequenceIn_] + columnLength[sequenceIn_]; i++) {
4208	int iRow = row[i];
4209	rowArray_[5]->quickAddNonZero(iRow, elementByColumn[i]scale rowScale_[iRow]*change);
4210	}
4211	}
4212	} else {
4213	rowArray_[5]->insert(sequenceIn_-numberColumns_,-change);
4214	}
4215	needed++;
4216	}
4217	if (getStatus(sequenceOut_)==atLowerBound)
4218	change=lowerChange[sequenceOut_];
4219	else
4220	change=upperChange[sequenceOut_];
4221	if (change) {
4222	if (sequenceOut_<numberColumns_) {
4223	if (!rowScale_) {
4224	for (CoinBigIndex i = columnStart[sequenceOut_];
4225	i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4226	rowArray_[5]->quickAddNonZero(row[i], elementByColumn[i]*change);
4227	}
4228	} else {
4229	// apply scaling
4230	double scale = columnScale_[sequenceOut_];
4231	for (CoinBigIndex i = columnStart[sequenceOut_];
4232	i < columnStart[sequenceOut_] + columnLength[sequenceOut_]; i++) {
4233	int iRow = row[i];
4234	rowArray_[5]->quickAddNonZero(iRow, elementByColumn[i]scale rowScale_[iRow]*change);
4235	}
4236	}
4237	} else {
4238	rowArray_[5]->quickAddNonZero(sequenceOut_-numberColumns_,-change);
4239	}
4240	needed++;
4241	}
4242	//printf("seqin %d seqout %d needed %d\n",
4243	// sequenceIn_,sequenceOut_,needed);
4244	if (needed) {
4245	// ftran it
4246	factorization_->updateColumn(rowArray_[0], rowArray_[5]);
4247	// add
4248	double * array5 = rowArray_[5]->denseVector();
4249	int * index5 = rowArray_[5]->getIndices();
4250	int number5 = rowArray_[5]->getNumElements();
4251	for (int i = 0; i < number5; i++) {
4252	int iPivot = index5[i];
4253	rowArray_[4]->quickAddNonZero(iPivot,array5[iPivot]);
4254	array5[iPivot]=0.0;
4255	}
4256	rowArray_[5]->setNumElements(0);
4257	}
4258	}
4259	pivotRow_ = -1;
4260	const int * index = rowArray_[4]->getIndices();
4261	int number = rowArray_[4]->getNumElements();
4262	int * lowerList2 = (reinterpret_cast<int >(lower_+4numberTotal))+2;
4263	char * markDone = reinterpret_cast<char *>(lowerList2+numberTotal);
4264	#if 0
4265	for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
4266	//int iPivot = index[iRow];
4267	iSequence = pivotVariable_[iPivot];
4268	// solution value will be sol - theta*alpha
4269	// bounds will be bounds + change *theta
4270	double currentSolution = solution_[iSequence];
4271	double currentLower = lower_[iSequence];
4272	double currentUpper = upper_[iSequence];
4273	double alpha = array[iPivot];
4274	assert (currentSolution >= currentLower - 100.0*primalTolerance_);
4275	assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
4276	double thetaCoefficient;
4277	double hitsLower = COIN_DBL_MAX;
4278	thetaCoefficient = lowerChange[iSequence] + alpha;
4279	if (thetaCoefficient > 1.0e-8)
4280	hitsLower = (currentSolution - currentLower) / thetaCoefficient;
4281	//if (hitsLower < 0.0) {
4282	// does not hit - but should we check further
4283	// hitsLower = COIN_DBL_MAX;
4284	//}
4285	double hitsUpper = COIN_DBL_MAX;
4286	thetaCoefficient = upperChange[iSequence] + alpha;
4287	if (thetaCoefficient < -1.0e-8)
4288	hitsUpper = (currentSolution - currentUpper) / thetaCoefficient;
4289	//if (hitsUpper < 0.0) {
4290	// does not hit - but should we check further
4291	// hitsUpper = COIN_DBL_MAX;
4292	//}
4293	if (CoinMin(hitsLower, hitsUpper) < theta_) {
4294	theta_ = CoinMin(hitsLower, hitsUpper);
4295	toLower = hitsLower < hitsUpper;
4296	pivotRow_ = iPivot;
4297	}
4298	}
4299	#else
4300	// first ones with alpha
4301	for (int i=0;i<number;i++) {
4302	int iPivot=index[i];
4303	iSequence = pivotVariable_[iPivot];
4304	assert(!markDone[iSequence]);
4305	markDone[iSequence]=1;
4306	// solution value will be sol - theta*alpha
4307	// bounds will be bounds + change *theta
4308	double currentSolution = solution_[iSequence];
4309	double currentLower = lower_[iSequence];
4310	double currentUpper = upper_[iSequence];
4311	double alpha = array[iPivot];
4312	assert (currentSolution >= currentLower - 100.0*primalTolerance_);
4313	assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
4314	double thetaCoefficient;
4315	thetaCoefficient = lowerChange[iSequence] + alpha;
4316	if (thetaCoefficient > 1.0e-8) {
4317	double gap=currentSolution-currentLower;
4318	if (thetaCoefficient*theta_>gap) {
4319	theta_ = gap/thetaCoefficient;
4320	toLower=true;
4321	pivotRow_=iPivot;
4322	}
4323	}
4324	thetaCoefficient = upperChange[iSequence] + alpha;
4325	if (thetaCoefficient < -1.0e-8) {
4326	double gap=currentSolution-currentUpper; //negative
4327	if (thetaCoefficient*theta_<gap) {
4328	theta_ = gap/thetaCoefficient;
4329	toLower=false;
4330	pivotRow_=iPivot;
4331	}
4332	}
4333	}
4334	// now others
4335	int nLook=lowerList[-1];
4336	for (int i=0;i<nLook;i++) {
4337	int iSequence = lowerList[i];
4338	if (getColumnStatus(iSequence)==basic&&!markDone[iSequence]) {
4339	double currentSolution = solution_[iSequence];
4340	double currentLower = lower_[iSequence];
4341	assert (currentSolution >= currentLower - 100.0*primalTolerance_);
4342	double thetaCoefficient = lowerChange[iSequence];
4343	if (thetaCoefficient > 0.0) {
4344	double gap=currentSolution-currentLower;
4345	if (thetaCoefficient*theta_>gap) {
4346	theta_ = gap/thetaCoefficient;
4347	toLower=true;
4348	pivotRow_ = -2-iSequence;
4349	}
4350	}
4351	}
4352	}
4353	nLook=upperList[-1];
4354	for (int i=0;i<nLook;i++) {
4355	int iSequence = upperList[i];
4356	if (getColumnStatus(iSequence)==basic&&!markDone[iSequence]) {
4357	double currentSolution = solution_[iSequence];
4358	double currentUpper = upper_[iSequence];
4359	assert (currentSolution <= currentUpper + 100.0*primalTolerance_);
4360	double thetaCoefficient = upperChange[iSequence];
4361	if (thetaCoefficient < 0) {
4362	double gap=currentSolution-currentUpper; //negative
4363	if (thetaCoefficient*theta_<gap) {
4364	theta_ = gap/thetaCoefficient;
4365	toLower=false;
4366	pivotRow_ = -2-iSequence;
4367	}
4368	}
4369	}
4370	}
4371	if (pivotRow_<-1) {
4372	// find
4373	int iSequence = -pivotRow_-2;
4374	for (int iPivot = 0; iPivot < numberRows_; iPivot++) {
4375	if (iSequence == pivotVariable_[iPivot]) {
4376	pivotRow_=iPivot;
4377	break;
4378	}
4379	}
4380	assert (pivotRow_>=0);
4381	}
4382	#endif
4383	theta_ = CoinMax(theta_,0.0);
4384	if (theta_>1.0e-15) {
4385	// update solution
4386	for (int iRow = 0; iRow < number; iRow++) {
4387	int iPivot = index[iRow];
4388	iSequence = pivotVariable_[iPivot];
4389	markDone[iSequence]=0;
4390	// solution value will be sol - theta*alpha
4391	double alpha = array[iPivot];
4392	solution_[iSequence] -= theta_ * alpha;
4393	}
4394	} else {
4395	for (int iRow = 0; iRow < number; iRow++) {
4396	int iPivot = index[iRow];
4397	iSequence = pivotVariable_[iPivot];
4398	markDone[iSequence]=0;
4399	}
4400	}
4401	#if 0
4402	for (int i=0;i<numberTotal;i++)
4403	assert(!markDone[i]);
4404	#endif
4405	if (pivotRow_ >= 0) {
4406	sequenceOut_ = pivotVariable_[pivotRow_];
4407	valueOut_ = solution_[sequenceOut_];
4408	lowerOut_ = lower_[sequenceOut_]+theta_*lowerChange[sequenceOut_];
4409	upperOut_ = upper_[sequenceOut_]+theta_*upperChange[sequenceOut_];
4410	if (!toLower) {
4411	directionOut_ = -1;
4412	dualOut_ = valueOut_ - upperOut_;
4413	} else {
4414	directionOut_ = 1;
4415	dualOut_ = lowerOut_ - valueOut_;
4416	}
4417	return 0;
4418	} else {
4419	//theta_=0.0;
4420	return -1;
4421	}
4422	}
4423	// Restores bound to original bound
4424	void
4425	ClpSimplexOther::originalBound(int iSequence, double theta,
4426	const double * lowerChange,
4427	const double * upperChange)
4428	{
4429	if (getFakeBound(iSequence) != noFake) {
4430	numberFake_--;
4431	setFakeBound(iSequence, noFake);
4432	if (iSequence >= numberColumns_) {
4433	// rows
4434	int iRow = iSequence - numberColumns_;
4435	rowLowerWork_[iRow] = rowLower_[iRow]+theta*lowerChange[iSequence];
4436	rowUpperWork_[iRow] = rowUpper_[iRow]+theta*upperChange[iSequence];
4437	if (rowScale_) {
4438	if (rowLowerWork_[iRow] > -1.0e50)
4439	rowLowerWork_[iRow] = rowScale_[iRow] rhsScale_;
4440	if (rowUpperWork_[iRow] < 1.0e50)
4441	rowUpperWork_[iRow] = rowScale_[iRow] rhsScale_;
4442	} else if (rhsScale_ != 1.0) {
4443	if (rowLowerWork_[iRow] > -1.0e50)
4444	rowLowerWork_[iRow] *= rhsScale_;
4445	if (rowUpperWork_[iRow] < 1.0e50)
4446	rowUpperWork_[iRow] *= rhsScale_;
4447	}
4448	} else {
4449	// columns
4450	columnLowerWork_[iSequence] = columnLower_[iSequence]+theta*lowerChange[iSequence];
4451	columnUpperWork_[iSequence] = columnUpper_[iSequence]+theta*upperChange[iSequence];
4452	if (rowScale_) {
4453	double multiplier = 1.0 * inverseColumnScale_[iSequence];
4454	if (columnLowerWork_[iSequence] > -1.0e50)
4455	columnLowerWork_[iSequence] = multiplier rhsScale_;
4456	if (columnUpperWork_[iSequence] < 1.0e50)
4457	columnUpperWork_[iSequence] = multiplier rhsScale_;
4458	} else if (rhsScale_ != 1.0) {
4459	if (columnLowerWork_[iSequence] > -1.0e50)
4460	columnLowerWork_[iSequence] *= rhsScale_;
4461	if (columnUpperWork_[iSequence] < 1.0e50)
4462	columnUpperWork_[iSequence] *= rhsScale_;
4463	}
4464	}
4465	}
4466	}
4467	/* Expands out all possible combinations for a knapsack
4468	If buildObj NULL then just computes space needed - returns number elements
4469	On entry numberOutput is maximum allowed, on exit it is number needed or
4470	-1 (as will be number elements) if maximum exceeded. numberOutput will have at
4471	least space to return values which reconstruct input.
4472	Rows returned will be original rows but no entries will be returned for
4473	any rows all of whose entries are in knapsack. So up to user to allow for this.
4474	If reConstruct >=0 then returns number of entrie which make up item "reConstruct"
4475	in expanded knapsack. Values in buildRow and buildElement;
4476	*/
4477	int
4478	ClpSimplexOther::expandKnapsack(int knapsackRow, int & numberOutput,
4479	double * buildObj, CoinBigIndex * buildStart,
4480	int * buildRow, double * buildElement, int reConstruct) const
4481	{
4482	int iRow;
4483	int iColumn;
4484	// Get column copy
4485	CoinPackedMatrix * columnCopy = matrix();
4486	// Get a row copy in standard format
4487	CoinPackedMatrix matrixByRow;
4488	matrixByRow.reverseOrderedCopyOf(*columnCopy);
4489	const double * elementByRow = matrixByRow.getElements();
4490	const int * column = matrixByRow.getIndices();
4491	const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
4492	const int * rowLength = matrixByRow.getVectorLengths();
4493	CoinBigIndex j;
4494	int * whichColumn = new int [numberColumns_];
4495	int * whichRow = new int [numberRows_];
4496	int numJ = 0;
4497	// Get what other columns can compensate for
4498	double * lo = new double [numberRows_];
4499	double * high = new double [numberRows_];
4500	{
4501	// Use to get tight column bounds
4502	ClpSimplex tempModel(*this);
4503	tempModel.tightenPrimalBounds(0.0, 0, true);
4504	// Now another model without knapsacks
4505	int nCol = 0;
4506	for (iRow = 0; iRow < numberRows_; iRow++) {
4507	whichRow[iRow] = iRow;
4508	}
4509	for (iColumn = 0; iColumn < numberColumns_; iColumn++)
4510	whichColumn[iColumn] = -1;
4511	for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
4512	int iColumn = column[j];
4513	if (columnUpper_[iColumn] > columnLower_[iColumn]) {
4514	whichColumn[iColumn] = 0;
4515	} else {
4516	assert (!columnLower_[iColumn]); // fix later
4517	}
4518	}
4519	for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
4520	if (whichColumn[iColumn] < 0)
4521	whichColumn[nCol++] = iColumn;
4522	}
4523	ClpSimplex tempModel2(&tempModel, numberRows_, whichRow, nCol, whichColumn, false, false, false);
4524	// Row copy
4525	CoinPackedMatrix matrixByRow;
4526	matrixByRow.reverseOrderedCopyOf(*tempModel2.matrix());
4527	const double * elementByRow = matrixByRow.getElements();
4528	const int * column = matrixByRow.getIndices();
4529	const CoinBigIndex * rowStart = matrixByRow.getVectorStarts();
4530	const int * rowLength = matrixByRow.getVectorLengths();
4531	const double * columnLower = tempModel2.getColLower();
4532	const double * columnUpper = tempModel2.getColUpper();
4533	for (iRow = 0; iRow < numberRows_; iRow++) {
4534	lo[iRow] = -COIN_DBL_MAX;
4535	high[iRow] = COIN_DBL_MAX;
4536	if (rowLower_[iRow] > -1.0e20 \|\| rowUpper_[iRow] < 1.0e20) {
4537
4538	// possible row
4539	int infiniteUpper = 0;
4540	int infiniteLower = 0;
4541	double maximumUp = 0.0;
4542	double maximumDown = 0.0;
4543	CoinBigIndex rStart = rowStart[iRow];
4544	CoinBigIndex rEnd = rowStart[iRow] + rowLength[iRow];
4545	CoinBigIndex j;
4546	// Compute possible lower and upper ranges
4547
4548	for (j = rStart; j < rEnd; ++j) {
4549	double value = elementByRow[j];
4550	iColumn = column[j];
4551	if (value > 0.0) {
4552	if (columnUpper[iColumn] >= 1.0e20) {
4553	++infiniteUpper;
4554	} else {
4555	maximumUp += columnUpper[iColumn] * value;
4556	}
4557	if (columnLower[iColumn] <= -1.0e20) {
4558	++infiniteLower;
4559	} else {
4560	maximumDown += columnLower[iColumn] * value;
4561	}
4562	} else if (value < 0.0) {
4563	if (columnUpper[iColumn] >= 1.0e20) {
4564	++infiniteLower;
4565	} else {
4566	maximumDown += columnUpper[iColumn] * value;
4567	}
4568	if (columnLower[iColumn] <= -1.0e20) {
4569	++infiniteUpper;
4570	} else {
4571	maximumUp += columnLower[iColumn] * value;
4572	}
4573	}
4574	}
4575	// Build in a margin of error
4576	maximumUp += 1.0e-8 * fabs(maximumUp) + 1.0e-7;
4577	maximumDown -= 1.0e-8 * fabs(maximumDown) + 1.0e-7;
4578	// we want to save effective rhs
4579	double up = (infiniteUpper) ? COIN_DBL_MAX : maximumUp;
4580	double down = (infiniteLower) ? -COIN_DBL_MAX : maximumDown;
4581	if (up == COIN_DBL_MAX \|\| rowLower_[iRow] == -COIN_DBL_MAX) {
4582	// However low we go it doesn't matter
4583	lo[iRow] = -COIN_DBL_MAX;
4584	} else {
4585	// If we go below this then can not be feasible
4586	lo[iRow] = rowLower_[iRow] - up;
4587	}
4588	if (down == -COIN_DBL_MAX \|\| rowUpper_[iRow] == COIN_DBL_MAX) {
4589	// However high we go it doesn't matter
4590	high[iRow] = COIN_DBL_MAX;
4591	} else {
4592	// If we go above this then can not be feasible
4593	high[iRow] = rowUpper_[iRow] - down;
4594	}
4595	}
4596	}
4597	}
4598	numJ = 0;
4599	for (iColumn = 0; iColumn < numberColumns_; iColumn++)
4600	whichColumn[iColumn] = -1;
4601	int * markRow = new int [numberRows_];
4602	for (iRow = 0; iRow < numberRows_; iRow++)
4603	markRow[iRow] = 1;
4604	for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
4605	int iColumn = column[j];
4606	if (columnUpper_[iColumn] > columnLower_[iColumn]) {
4607	whichColumn[iColumn] = numJ;
4608	numJ++;
4609	}
4610	}
4611	/* mark rows
4612	-n in knapsack and n other variables
4613	1 no entries
4614	n+1000 not involved in knapsack but n entries
4615	0 only in knapsack
4616	*/
4617	for (iRow = 0; iRow < numberRows_; iRow++) {
4618	int type = 1;
4619	for (j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
4620	int iColumn = column[j];
4621	if (whichColumn[iColumn] >= 0) {
4622	if (type == 1) {
4623	type = 0;
4624	} else if (type > 0) {
4625	assert (type > 1000);
4626	type = -(type - 1000);
4627	}
4628	} else if (type == 1) {
4629	type = 1001;
4630	} else if (type < 0) {
4631	type --;
4632	} else if (type == 0) {
4633	type = -1;
4634	} else {
4635	assert (type > 1000);
4636	type++;
4637	}
4638	}
4639	markRow[iRow] = type;
4640	if (type < 0 && type > -30 && false)
4641	printf("markrow on row %d is %d\n", iRow, markRow[iRow]);
4642	}
4643	int * bound = new int [numberColumns_+1];
4644	int * stack = new int [numberColumns_+1];
4645	int * flip = new int [numberColumns_+1];
4646	double * offset = new double[numberColumns_+1];
4647	double * size = new double [numberColumns_+1];
4648	double * rhsOffset = new double[numberRows_];
4649	int * build = new int[numberColumns_];
4650	int maxNumber = numberOutput;
4651	numJ = 0;
4652	double minSize = rowLower_[knapsackRow];
4653	double maxSize = rowUpper_[knapsackRow];
4654	double knapsackOffset = 0.0;
4655	for (j = rowStart[knapsackRow]; j < rowStart[knapsackRow] + rowLength[knapsackRow]; j++) {
4656	int iColumn = column[j];
4657	double lowerColumn = columnLower_[iColumn];
4658	double upperColumn = columnUpper_[iColumn];
4659	if (lowerColumn == upperColumn)
4660	continue;
4661	double gap = upperColumn - lowerColumn;
4662	if (gap > 1.0e8)
4663	gap = 1.0e8;
4664	assert (fabs(floor(gap + 0.5) - gap) < 1.0e-5);
4665	whichColumn[numJ] = iColumn;
4666	bound[numJ] = static_cast<int> (gap);
4667	if (elementByRow[j] > 0.0) {
4668	flip[numJ] = 1;
4669	offset[numJ] = lowerColumn;
4670	size[numJ++] = elementByRow[j];
4671	} else {
4672	flip[numJ] = -1;
4673	offset[numJ] = upperColumn;
4674	size[numJ++] = -elementByRow[j];
4675	lowerColumn = upperColumn;
4676	}
4677	knapsackOffset += elementByRow[j] * lowerColumn;
4678	}
4679	int jRow;
4680	for (iRow = 0; iRow < numberRows_; iRow++)
4681	whichRow[iRow] = iRow;
4682	ClpSimplex smallModel(this, numberRows_, whichRow, numJ, whichColumn, true, true, true);
4683	// modify rhs to allow for nonzero lower bounds
4684	//double * rowLower = smallModel.rowLower();
4685	//double * rowUpper = smallModel.rowUpper();
4686	//const double * columnLower = smallModel.columnLower();
4687	//const double * columnUpper = smallModel.columnUpper();
4688	const CoinPackedMatrix * matrix = smallModel.matrix();
4689	const double * element = matrix->getElements();
4690	const int * row = matrix->getIndices();
4691	const CoinBigIndex * columnStart = matrix->getVectorStarts();
4692	const int * columnLength = matrix->getVectorLengths();
4693	const double * objective = smallModel.objective();
4694	//double objectiveOffset=0.0;
4695	// would use for fixed?
4696	CoinZeroN(rhsOffset, numberRows_);
4697	double * rowActivity = smallModel.primalRowSolution();
4698	CoinZeroN(rowActivity, numberRows_);
4699	maxSize -= knapsackOffset;
4700	minSize -= knapsackOffset;
4701	// now generate
4702	int i;
4703	int iStack = numJ;
4704	for (i = 0; i < numJ; i++) {
4705	stack[i] = 0;
4706	}
4707	double tooMuch = 10.0 * maxSize + 10000;
4708	stack[numJ] = 1;
4709	size[numJ] = tooMuch;
4710	bound[numJ] = 0;
4711	double sum = tooMuch;
4712	// allow for all zero being OK
4713	stack[numJ-1] = -1;
4714	sum -= size[numJ-1];
4715	numberOutput = 0;
4716	int nelCreate = 0;
4717	/* typeRun is - 0 for initial sizes
4718	1 for build
4719	2 for reconstruct
4720	*/
4721	int typeRun = buildObj ? 1 : 0;
4722	if (reConstruct >= 0) {
4723	assert (buildRow && buildElement);
4724	typeRun = 2;
4725	}
4726	if (typeRun == 1)
4727	buildStart[0] = 0;
4728	while (iStack >= 0) {
4729	if (sum >= minSize && sum <= maxSize) {
4730	double checkSize = 0.0;
4731	bool good = true;
4732	int nRow = 0;
4733	double obj = 0.0;
4734	CoinZeroN(rowActivity, numberRows_);
4735	for (iColumn = 0; iColumn < numJ; iColumn++) {
4736	int iValue = stack[iColumn];
4737	if (iValue > bound[iColumn]) {
4738	good = false;
4739	break;
4740	} else {
4741	double realValue = offset[iColumn] + flip[iColumn] * iValue;
4742	if (realValue) {
4743	obj += objective[iColumn] * realValue;
4744	for (CoinBigIndex j = columnStart[iColumn];
4745	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4746	double value = element[j] * realValue;
4747	int kRow = row[j];
4748	if (rowActivity[kRow]) {
4749	rowActivity[kRow] += value;
4750	if (!rowActivity[kRow])
4751	rowActivity[kRow] = 1.0e-100;
4752	} else {
4753	build[nRow++] = kRow;
4754	rowActivity[kRow] = value;
4755	}
4756	}
4757	}
4758	}
4759	}
4760	if (good) {
4761	for (jRow = 0; jRow < nRow; jRow++) {
4762	int kRow = build[jRow];
4763	double value = rowActivity[kRow];
4764	if (value > high[kRow] \|\| value < lo[kRow]) {
4765	good = false;
4766	break;
4767	}
4768	}
4769	}
4770	if (good) {
4771	if (typeRun == 1) {
4772	buildObj[numberOutput] = obj;
4773	for (jRow = 0; jRow < nRow; jRow++) {
4774	int kRow = build[jRow];
4775	double value = rowActivity[kRow];
4776	if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
4777	buildElement[nelCreate] = value;
4778	buildRow[nelCreate++] = kRow;
4779	}
4780	}
4781	buildStart[numberOutput+1] = nelCreate;
4782	} else if (!typeRun) {
4783	for (jRow = 0; jRow < nRow; jRow++) {
4784	int kRow = build[jRow];
4785	double value = rowActivity[kRow];
4786	if (markRow[kRow] < 0 && fabs(value) > 1.0e-13) {
4787	nelCreate++;
4788	}
4789	}
4790	}
4791	if (typeRun == 2 && reConstruct == numberOutput) {
4792	// build and exit
4793	nelCreate = 0;
4794	for (iColumn = 0; iColumn < numJ; iColumn++) {
4795	int iValue = stack[iColumn];
4796	double realValue = offset[iColumn] + flip[iColumn] * iValue;
4797	if (realValue) {
4798	buildRow[nelCreate] = whichColumn[iColumn];
4799	buildElement[nelCreate++] = realValue;
4800	}
4801	}
4802	numberOutput = 1;
4803	for (i = 0; i < numJ; i++) {
4804	bound[i] = 0;
4805	}
4806	break;
4807	}
4808	numberOutput++;
4809	if (numberOutput > maxNumber) {
4810	nelCreate = -numberOutput;
4811	numberOutput = -1;
4812	for (i = 0; i < numJ; i++) {
4813	bound[i] = 0;
4814	}
4815	break;
4816	} else if (typeRun == 1 && numberOutput == maxNumber) {
4817	// On second run
4818	for (i = 0; i < numJ; i++) {
4819	bound[i] = 0;
4820	}
4821	break;
4822	}
4823	for (int j = 0; j < numJ; j++) {
4824	checkSize += stack[j] * size[j];
4825	}
4826	assert (fabs(sum - checkSize) < 1.0e-3);
4827	}
4828	for (jRow = 0; jRow < nRow; jRow++) {
4829	int kRow = build[jRow];
4830	rowActivity[kRow] = 0.0;
4831	}
4832	}
4833	if (sum > maxSize \|\| stack[iStack] > bound[iStack]) {
4834	sum -= size[iStack] * stack[iStack];
4835	stack[iStack--] = 0;
4836	if (iStack >= 0) {
4837	stack[iStack] ++;
4838	sum += size[iStack];
4839	}
4840	} else {
4841	// must be less
4842	// add to last possible
4843	iStack = numJ - 1;
4844	sum += size[iStack];
4845	stack[iStack]++;
4846	}
4847	}
4848	//printf("%d will be created\n",numberOutput);
4849	delete [] whichColumn;
4850	delete [] whichRow;
4851	delete [] bound;
4852	delete [] stack;
4853	delete [] flip;
4854	delete [] size;
4855	delete [] offset;
4856	delete [] rhsOffset;
4857	delete [] build;
4858	delete [] markRow;
4859	delete [] lo;
4860	delete [] high;
4861	return nelCreate;
4862	}
4863	// Quick try at cleaning up duals if postsolve gets wrong
4864	void
4865	ClpSimplexOther::cleanupAfterPostsolve()
4866	{
4867	// First mark singleton equality rows
4868	char * mark = new char [ numberRows_];
4869	memset(mark, 0, numberRows_);
4870	const int * row = matrix_->getIndices();
4871	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4872	const int * columnLength = matrix_->getVectorLengths();
4873	const double * element = matrix_->getElements();
4874	for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
4875	for (CoinBigIndex j = columnStart[iColumn];
4876	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4877	int iRow = row[j];
4878	if (mark[iRow])
4879	mark[iRow] = 2;
4880	else
4881	mark[iRow] = 1;
4882	}
4883	}
4884	// for now just == rows
4885	for (int iRow = 0; iRow < numberRows_; iRow++) {
4886	if (rowUpper_[iRow] > rowLower_[iRow])
4887	mark[iRow] = 3;
4888	}
4889	double dualTolerance = dblParam_[ClpDualTolerance];
4890	double primalTolerance = dblParam_[ClpPrimalTolerance];
4891	int numberCleaned = 0;
4892	double maxmin = optimizationDirection_;
4893	for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
4894	double dualValue = reducedCost_[iColumn] * maxmin;
4895	double primalValue = columnActivity_[iColumn];
4896	double lower = columnLower_[iColumn];
4897	double upper = columnUpper_[iColumn];
4898	int way = 0;
4899	switch(getColumnStatus(iColumn)) {
4900
4901	case basic:
4902	// dual should be zero
4903	if (dualValue > dualTolerance) {
4904	way = -1;
4905	} else if (dualValue < -dualTolerance) {
4906	way = 1;
4907	}
4908	break;
4909	case ClpSimplex::isFixed:
4910	break;
4911	case atUpperBound:
4912	// dual should not be positive
4913	if (dualValue > dualTolerance) {
4914	way = -1;
4915	}
4916	break;
4917	case atLowerBound:
4918	// dual should not be negative
4919	if (dualValue < -dualTolerance) {
4920	way = 1;
4921	}
4922	break;
4923	case superBasic:
4924	case isFree:
4925	if (primalValue < upper - primalTolerance) {
4926	// dual should not be negative
4927	if (dualValue < -dualTolerance) {
4928	way = 1;
4929	}
4930	}
4931	if (primalValue > lower + primalTolerance) {
4932	// dual should not be positive
4933	if (dualValue > dualTolerance) {
4934	way = -1;
4935	}
4936	}
4937	break;
4938	}
4939	if (way) {
4940	// see if can find singleton row
4941	for (CoinBigIndex j = columnStart[iColumn];
4942	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4943	int iRow = row[j];
4944	if (mark[iRow] == 1) {
4945	double value = element[j];
4946	// dj - addDual*value == 0.0
4947	double addDual = dualValue / value;
4948	dual_[iRow] += addDual;
4949	reducedCost_[iColumn] = 0.0;
4950	numberCleaned++;
4951	break;
4952	}
4953	}
4954	}
4955	}
4956	delete [] mark;
4957	#ifdef CLP_INVESTIGATE
4958	printf("cleanupAfterPostsolve cleaned up %d columns\n", numberCleaned);
4959	#endif
4960	// Redo
4961	memcpy(reducedCost_, this->objective(), numberColumns_ * sizeof(double));
4962	matrix_->transposeTimes(-1.0, dual_, reducedCost_);
4963	checkSolutionInternal();
4964	}
4965	// Returns gub version of model or NULL
4966	ClpSimplex *
4967	ClpSimplexOther::gubVersion(int * whichRows, int * whichColumns,
4968	int neededGub,
4969	int factorizationFrequency)
4970	{
4971	// find gub
4972	int numberRows = this->numberRows();
4973	int numberColumns = this->numberColumns();
4974	int iRow, iColumn;
4975	int * columnIsGub = new int [numberColumns];
4976	const double * columnLower = this->columnLower();
4977	const double * columnUpper = this->columnUpper();
4978	int numberFixed=0;
4979	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
4980	if (columnUpper[iColumn] == columnLower[iColumn]) {
4981	columnIsGub[iColumn]=-2;
4982	numberFixed++;
4983	} else if (columnLower[iColumn]>=0) {
4984	columnIsGub[iColumn]=-1;
4985	} else {
4986	columnIsGub[iColumn]=-3;
4987	}
4988	}
4989	CoinPackedMatrix * matrix = this->matrix();
4990	// get row copy
4991	CoinPackedMatrix rowCopy = *matrix;
4992	rowCopy.reverseOrdering();
4993	const int * column = rowCopy.getIndices();
4994	const int * rowLength = rowCopy.getVectorLengths();
4995	const CoinBigIndex * rowStart = rowCopy.getVectorStarts();
4996	const double * element = rowCopy.getElements();
4997	int numberNonGub = 0;
4998	int numberEmpty = numberRows;
4999	int * rowIsGub = new int [numberRows];
5000	int smallestGubRow=-1;
5001	int count=numberColumns+1;
5002	double * rowLower = this->rowLower();
5003	double * rowUpper = this->rowUpper();
5004	// make sure we can get rid of upper bounds
5005	double * fixedRow = new double [numberRows];
5006	for (iRow = 0 ; iRow < numberRows ; iRow++) {
5007	double sumFixed=0.0;
5008	for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
5009	int iColumn = column[j];
5010	double value = columnLower[iColumn];
5011	if (value)
5012	sumFixed += element[j] * value;
5013	}
5014	fixedRow[iRow]=rowUpper[iRow]-sumFixed;
5015	}
5016	for (iRow = numberRows - 1; iRow >= 0; iRow--) {
5017	bool gubRow = true;
5018	int numberInRow=0;
5019	double sumFixed=0.0;
5020	double gap = fixedRow[iRow]-1.0e-12;
5021	for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
5022	int iColumn = column[j];
5023	if (columnIsGub[iColumn]!=-2) {
5024	if (element[j] != 1.0\|\|columnIsGub[iColumn]==-3\|\|
5025	columnUpper[iColumn]-columnLower[iColumn]<gap) {
5026	gubRow = false;
5027	break;
5028	} else {
5029	numberInRow++;
5030	if (columnIsGub[iColumn] >= 0) {
5031	gubRow = false;
5032	break;
5033	}
5034	}
5035	} else {
5036	sumFixed += columnLower[iColumn]*element[j];
5037	}
5038	}
5039	if (!gubRow) {
5040	whichRows[numberNonGub++] = iRow;
5041	rowIsGub[iRow] = -1;
5042	} else if (numberInRow) {
5043	if (numberInRow<count) {
5044	count = numberInRow;
5045	smallestGubRow=iRow;
5046	}
5047	for (int j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
5048	int iColumn = column[j];
5049	if (columnIsGub[iColumn]!=-2)
5050	columnIsGub[iColumn] = iRow;
5051	}
5052	rowIsGub[iRow] = 0;
5053	} else {
5054	// empty row!
5055	whichRows[--numberEmpty] = iRow;
5056	rowIsGub[iRow] = -2;
5057	if (sumFixed>rowUpper[iRow]+1.0e-4\|\|
5058	sumFixed<rowLower[iRow]-1.0e-4) {
5059	fprintf(stderr,"******** No infeasible empty rows - please!\n");
5060	abort();
5061	}
5062	}
5063	}
5064	delete [] fixedRow;
5065	char message[100];
5066	int numberGub = numberEmpty - numberNonGub;
5067	if (numberGub >= neededGub) {
5068	sprintf(message,"%d gub rows", numberGub);
5069	handler_->message(CLP_GENERAL2, messages_)
5070	<< message << CoinMessageEol;
5071	int numberNormal = 0;
5072	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5073	if (columnIsGub[iColumn] < 0 && columnIsGub[iColumn] !=-2) {
5074	whichColumns[numberNormal++] = iColumn;
5075	}
5076	}
5077	if (!numberNormal) {
5078	sprintf(message,"Putting back one gub row to make non-empty");
5079	handler_->message(CLP_GENERAL2, messages_)
5080	<< message << CoinMessageEol;
5081	rowIsGub[smallestGubRow]=-1;
5082	whichRows[numberNonGub++] = smallestGubRow;
5083	for (int j = rowStart[smallestGubRow];
5084	j < rowStart[smallestGubRow] + rowLength[smallestGubRow]; j++) {
5085	int iColumn = column[j];
5086	if (columnIsGub[iColumn]>=0) {
5087	columnIsGub[iColumn]=-4;
5088	whichColumns[numberNormal++] = iColumn;
5089	}
5090	}
5091	}
5092	std::sort(whichRows,whichRows+numberNonGub);
5093	std::sort(whichColumns,whichColumns+numberNormal);
5094	double * lower = CoinCopyOfArray(this->rowLower(),numberRows);
5095	double * upper = CoinCopyOfArray(this->rowUpper(),numberRows);
5096	// leave empty rows at end
5097	numberEmpty = numberRows-numberEmpty;
5098	const int * row = matrix->getIndices();
5099	const int * columnLength = matrix->getVectorLengths();
5100	const CoinBigIndex * columnStart = matrix->getVectorStarts();
5101	const double * elementByColumn = matrix->getElements();
5102	// Fixed at end
5103	int put2 = numberColumns-numberFixed;
5104	for (iColumn = 0; iColumn < numberColumns; iColumn++) {
5105	if (columnIsGub[iColumn] ==-2) {
5106	whichColumns[put2++] = iColumn;
5107	double value = columnLower[iColumn];
5108	for (int j = columnStart[iColumn];
5109	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
5110	int iRow = row[j];
5111	if (lower[iRow]>-1.0e20)
5112	lower[iRow] -= value*element[j];
5113	if (upper[iRow]<1.0e20)
5114	upper[iRow] -= value*element[j];
5115	}
5116	}
5117	}
5118	int put = numberNormal;
5119	ClpSimplex * model2 =
5120	new ClpSimplex(this, numberNonGub, whichRows , numberNormal, whichColumns);
5121	// scale
5122	double * scaleArray = new double [numberRows];
5123	for (int i=0;i<numberRows;i++) {
5124	scaleArray[i]=1.0;
5125	if (rowIsGub[i]==-1) {
5126	double largest = 1.0e-30;
5127	double smallest = 1.0e30;
5128	for (int j = rowStart[i]; j < rowStart[i] + rowLength[i]; j++) {
5129	int iColumn = column[j];
5130	if (columnIsGub[iColumn]!=-2) {
5131	double value =fabs(element[j]);
5132	largest = CoinMax(value,largest);
5133	smallest = CoinMin(value,smallest);
5134	}
5135	}
5136	double scale = CoinMax(0.001,1.0/sqrt(largest*smallest));
5137	scaleArray[i]=scale;
5138	if (lower[i]>-1.0e30)
5139	lower[i] *= scale;
5140	if (upper[i]<1.0e30)
5141	upper[i] *= scale;
5142	}
5143	}
5144	// scale partial matrix
5145	{
5146	CoinPackedMatrix * matrix = model2->matrix();
5147	const int * row = matrix->getIndices();
5148	const int * columnLength