Context Navigation

ClpSimplexOther.cpp @ 1860

Last change on this file since 1860 was 1860, checked in by forrest, 8 years ago
make test on maximum theta in parametrics more robust
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 354.8 KB

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