Context Navigation

ClpSimplexOther.cpp @ 2470

Last change on this file since 2470 was 2388, checked in by unxusr, 11 months ago
small change in parameter in guess
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 364.7 KB

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