Context Navigation

ClpSimplexOther.cpp @ 2128

Last change on this file since 2128 was 2128, checked in by forrest, 7 years ago
out mac warnings
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 363.0 KB

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