Context Navigation

source: trunk/Clp/src/ClpDynamicMatrix.cpp @ 754

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago
first version
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 73.7 KB

Line
1	// Copyright (C) 2004, International Business Machines
2	// Corporation and others. All Rights Reserved.
3
4
5	#include <cstdio>
6
7	#include "CoinPragma.hpp"
8	#include "CoinIndexedVector.hpp"
9	#include "CoinHelperFunctions.hpp"
10
11	#include "ClpSimplex.hpp"
12	#include "ClpFactorization.hpp"
13	#include "ClpQuadraticObjective.hpp"
14	#include "ClpNonLinearCost.hpp"
15	// at end to get min/max!
16	#include "ClpDynamicMatrix.hpp"
17	#include "ClpMessage.hpp"
18	//#define CLP_DEBUG
19	//#define CLP_DEBUG_PRINT
20	//#############################################################################
21	// Constructors / Destructor / Assignment
22	//#############################################################################
23
24	//-------------------------------------------------------------------
25	// Default Constructor
26	//-------------------------------------------------------------------
27	ClpDynamicMatrix::ClpDynamicMatrix ()
28	: ClpPackedMatrix(),
29	sumDualInfeasibilities_(0.0),
30	sumPrimalInfeasibilities_(0.0),
31	sumOfRelaxedDualInfeasibilities_(0.0),
32	sumOfRelaxedPrimalInfeasibilities_(0.0),
33	savedBestGubDual_(0.0),
34	savedBestSet_(0),
35	backToPivotRow_(NULL),
36	keyVariable_(NULL),
37	toIndex_(NULL),
38	fromIndex_(NULL),
39	numberSets_(0),
40	numberActiveSets_(0),
41	objectiveOffset_(0.0),
42	lowerSet_(NULL),
43	upperSet_(NULL),
44	status_(NULL),
45	model_(NULL),
46	firstAvailable_(0),
47	firstAvailableBefore_(0),
48	firstDynamic_(0),
49	lastDynamic_(0),
50	numberStaticRows_(0),
51	numberElements_(0),
52	numberDualInfeasibilities_(0),
53	numberPrimalInfeasibilities_(0),
54	noCheck_(-1),
55	infeasibilityWeight_(0.0),
56	numberGubColumns_(0),
57	maximumGubColumns_(0),
58	maximumElements_(0),
59	startSet_(NULL),
60	next_(NULL),
61	startColumn_(NULL),
62	row_(NULL),
63	element_(NULL),
64	cost_(NULL),
65	id_(NULL),
66	dynamicStatus_(NULL),
67	columnLower_(NULL),
68	columnUpper_(NULL)
69	{
70	setType(15);
71	}
72
73	//-------------------------------------------------------------------
74	// Copy constructor
75	//-------------------------------------------------------------------
76	ClpDynamicMatrix::ClpDynamicMatrix (const ClpDynamicMatrix & rhs)
77	: ClpPackedMatrix(rhs)
78	{
79	objectiveOffset_ = rhs.objectiveOffset_;
80	numberSets_ = rhs.numberSets_;
81	numberActiveSets_ = rhs.numberActiveSets_;
82	firstAvailable_ = rhs.firstAvailable_;
83	firstAvailableBefore_ = rhs.firstAvailableBefore_;
84	firstDynamic_ = rhs.firstDynamic_;
85	lastDynamic_ = rhs.lastDynamic_;
86	numberStaticRows_ = rhs.numberStaticRows_;
87	numberElements_ = rhs.numberElements_;
88	backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,lastDynamic_);
89	keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_);
90	toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_);
91	fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1-numberStaticRows_);
92	lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
93	upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
94	status_ = ClpCopyOfArray(rhs.status_,numberSets_);
95	model_ = rhs.model_;
96	sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
97	sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
98	sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
99	sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
100	numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
101	numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
102	savedBestGubDual_=rhs.savedBestGubDual_;
103	savedBestSet_=rhs.savedBestSet_;
104	noCheck_ = rhs.noCheck_;
105	infeasibilityWeight_=rhs.infeasibilityWeight_;
106	// Now secondary data
107	numberGubColumns_ = rhs.numberGubColumns_;
108	maximumGubColumns_=rhs.maximumGubColumns_;
109	maximumElements_ = rhs.maximumElements_;
110	startSet_ = ClpCopyOfArray(rhs.startSet_,numberSets_);
111	next_ = ClpCopyOfArray(rhs.next_,maximumGubColumns_);
112	startColumn_ = ClpCopyOfArray(rhs.startColumn_,maximumGubColumns_+1);
113	row_ = ClpCopyOfArray(rhs.row_,maximumElements_);
114	element_ = ClpCopyOfArray(rhs.element_,maximumElements_);
115	cost_ = ClpCopyOfArray(rhs.cost_,maximumGubColumns_);
116	id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_);
117	columnLower_ = ClpCopyOfArray(rhs.columnLower_,maximumGubColumns_);
118	columnUpper_ = ClpCopyOfArray(rhs.columnUpper_,maximumGubColumns_);
119	dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,maximumGubColumns_);
120	}
121
122	/* This is the real constructor*/
123	ClpDynamicMatrix::ClpDynamicMatrix(ClpSimplex * model, int numberSets,
124	int numberGubColumns, const int * starts,
125	const double * lower, const double * upper,
126	const CoinBigIndex * startColumn, const int * row,
127	const double * element, const double * cost,
128	const double * columnLower, const double * columnUpper,
129	const unsigned char * status,
130	const unsigned char * dynamicStatus)
131	: ClpPackedMatrix()
132	{
133	setType(15);
134	objectiveOffset_ = model->objectiveOffset();
135	model_ = model;
136	numberSets_ = numberSets;
137	numberGubColumns_ = numberGubColumns;
138	maximumGubColumns_=numberGubColumns_;
139	if (numberGubColumns_)
140	maximumElements_ = startColumn[numberGubColumns_];
141	else
142	maximumElements_ = 0;
143	startSet_ = new int [numberSets_];
144	next_ = new int [maximumGubColumns_];
145	// fill in startSet and next
146	int iSet;
147	if (numberGubColumns_) {
148	for (iSet=0;iSet<numberSets_;iSet++) {
149	int first = starts[iSet];
150	int last = starts[iSet+1]-1;
151	startSet_[iSet]=first;
152	for (int i=first;i<last;i++)
153	next_[i]=i+1;
154	next_[last]=-iSet-1;
155	}
156	}
157	int numberColumns = model->numberColumns();
158	int numberRows = model->numberRows();
159	numberStaticRows_ = numberRows;
160	savedBestGubDual_=0.0;
161	savedBestSet_=0;
162	// Number of columns needed
163	int frequency = model->factorizationFrequency();
164	int numberGubInSmall = numberRows + frequency + CoinMin(frequency,numberSets_)+4;
165	// for small problems this could be too big
166	//numberGubInSmall = CoinMin(numberGubInSmall,numberGubColumns_);
167	int numberNeeded = numberGubInSmall + numberColumns;
168	firstAvailable_ = numberColumns;
169	firstAvailableBefore_ = firstAvailable_;
170	firstDynamic_ = numberColumns;
171	lastDynamic_ = numberNeeded;
172	startColumn_ = ClpCopyOfArray(startColumn,numberGubColumns_+1);
173	if (!numberGubColumns_) {
174	startColumn_ = new CoinBigIndex [1];
175	startColumn_[0]=0;
176	}
177	CoinBigIndex numberElements = startColumn_[numberGubColumns_];
178	row_ = ClpCopyOfArray(row,numberElements);
179	element_ = new float[numberElements];
180	CoinBigIndex i;
181	for (i=0;i<numberElements;i++)
182	element_[i]=element[i];
183	cost_ = new float[numberGubColumns_];
184	for (i=0;i<numberGubColumns_;i++) {
185	cost_[i]=cost[i];
186	// I don't think I need sorted but ...
187	CoinSort_2(row_+startColumn_[i],row_+startColumn_[i+1],element_+startColumn_[i]);
188	}
189	if (columnLower) {
190	columnLower_ = new float[numberGubColumns_];
191	for (i=0;i<numberGubColumns_;i++)
192	columnLower_[i]=columnLower[i];
193	} else {
194	columnLower_=NULL;
195	}
196	if (columnUpper) {
197	columnUpper_ = new float[numberGubColumns_];
198	for (i=0;i<numberGubColumns_;i++)
199	columnUpper_[i]=columnUpper[i];
200	} else {
201	columnUpper_=NULL;
202	}
203	lowerSet_ = new float[numberSets_];
204	for (i=0;i<numberSets_;i++) {
205	if (lower[i]>-1.0e20)
206	lowerSet_[i]=lower[i];
207	else
208	lowerSet_[i] = -1.0e30;
209	}
210	upperSet_ = new float[numberSets_];
211	for (i=0;i<numberSets_;i++) {
212	if (upper[i]<1.0e20)
213	upperSet_[i]=upper[i];
214	else
215	upperSet_[i] = 1.0e30;
216	}
217	id_ = new int[numberGubInSmall];
218	for (i=0;i<numberGubInSmall;i++)
219	id_[i]=-1;
220	ClpPackedMatrix* originalMatrixA =
221	dynamic_cast< ClpPackedMatrix*>(model->clpMatrix());
222	assert (originalMatrixA);
223	CoinPackedMatrix * originalMatrix = originalMatrixA->getPackedMatrix();
224	originalMatrixA->setMatrixNull(); // so can be deleted safely
225	// guess how much space needed
226	double guess = numberElements;
227	guess /= (double) numberColumns;
228	guess = 2numberGubInSmall;
229	numberElements_ = (int) guess;
230	numberElements_ = CoinMin(numberElements_,numberElements)+originalMatrix->getNumElements();
231	matrix_ = originalMatrix;
232	zeroElements_ = false;
233	// resize model (matrix stays same)
234	int newRowSize = numberRows+CoinMin(numberSets_,CoinMax(frequency,numberRows))+1;
235	model->resize(newRowSize,numberNeeded);
236	for (i=numberRows;i<newRowSize;i++)
237	model->setRowStatus(i,ClpSimplex::basic);
238	if (columnUpper_) {
239	// set all upper bounds so we have enough space
240	double * columnUpper = model->columnUpper();
241	for(i=firstDynamic_;i<lastDynamic_;i++)
242	columnUpper[i]=1.0e10;
243	}
244	// resize matrix
245	// extra 1 is so can keep number of elements handy
246	originalMatrix->reserve(numberNeeded,numberElements_,true);
247	originalMatrix->reserve(numberNeeded+1,numberElements_,false);
248	originalMatrix->getMutableVectorStarts()[numberColumns]=originalMatrix->getNumElements();
249	originalMatrix->setDimensions(newRowSize,-1);
250	numberActiveColumns_=firstDynamic_;
251	// redo number of columns
252	numberColumns = matrix_->getNumCols();
253	backToPivotRow_ = new int[numberNeeded];
254	keyVariable_ = new int[numberSets_];
255	if (status) {
256	status_ = ClpCopyOfArray(status,numberSets_);
257	assert (dynamicStatus);
258	dynamicStatus_ = ClpCopyOfArray(dynamicStatus,numberGubColumns_);
259	} else {
260	status_= new unsigned char [numberSets_];
261	memset(status_,0,numberSets_);
262	int i;
263	for (i=0;i<numberSets_;i++) {
264	// make slack key
265	setStatus(i,ClpSimplex::basic);
266	}
267	dynamicStatus_ = new unsigned char [numberGubColumns_];
268	memset(dynamicStatus_,0,numberGubColumns_); // for clarity
269	for (i=0;i<numberGubColumns_;i++)
270	setDynamicStatus(i,atLowerBound);
271	}
272	toIndex_ = new int[numberSets_];
273	for (iSet=0;iSet<numberSets_;iSet++)
274	toIndex_[iSet]=-1;
275	fromIndex_ = new int [newRowSize-numberStaticRows_+1];
276	numberActiveSets_=0;
277	rhsOffset_=NULL;
278	if (numberGubColumns_) {
279	if (!status) {
280	gubCrash();
281	} else {
282	initialProblem();
283	}
284	}
285	noCheck_ = -1;
286	infeasibilityWeight_=0.0;
287	}
288
289	//-------------------------------------------------------------------
290	// Destructor
291	//-------------------------------------------------------------------
292	ClpDynamicMatrix::~ClpDynamicMatrix ()
293	{
294	delete [] backToPivotRow_;
295	delete [] keyVariable_;
296	delete [] toIndex_;
297	delete [] fromIndex_;
298	delete [] lowerSet_;
299	delete [] upperSet_;
300	delete [] status_;
301	delete [] startSet_;
302	delete [] next_;
303	delete [] startColumn_;
304	delete [] row_;
305	delete [] element_;
306	delete [] cost_;
307	delete [] id_;
308	delete [] dynamicStatus_;
309	delete [] columnLower_;
310	delete [] columnUpper_;
311	}
312
313	//----------------------------------------------------------------
314	// Assignment operator
315	//-------------------------------------------------------------------
316	ClpDynamicMatrix &
317	ClpDynamicMatrix::operator=(const ClpDynamicMatrix& rhs)
318	{
319	if (this != &rhs) {
320	ClpPackedMatrix::operator=(rhs);
321	delete [] backToPivotRow_;
322	delete [] keyVariable_;
323	delete [] toIndex_;
324	delete [] fromIndex_;
325	delete [] lowerSet_;
326	delete [] upperSet_;
327	delete [] status_;
328	delete [] startSet_;
329	delete [] next_;
330	delete [] startColumn_;
331	delete [] row_;
332	delete [] element_;
333	delete [] cost_;
334	delete [] id_;
335	delete [] dynamicStatus_;
336	delete [] columnLower_;
337	delete [] columnUpper_;
338	objectiveOffset_ = rhs.objectiveOffset_;
339	numberSets_ = rhs.numberSets_;
340	numberActiveSets_ = rhs.numberActiveSets_;
341	firstAvailable_ = rhs.firstAvailable_;
342	firstAvailableBefore_ = rhs.firstAvailableBefore_;
343	firstDynamic_ = rhs.firstDynamic_;
344	lastDynamic_ = rhs.lastDynamic_;
345	numberStaticRows_ = rhs.numberStaticRows_;
346	numberElements_ = rhs.numberElements_;
347	backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_,lastDynamic_);
348	keyVariable_ = ClpCopyOfArray(rhs.keyVariable_,numberSets_);
349	toIndex_ = ClpCopyOfArray(rhs.toIndex_,numberSets_);
350	fromIndex_ = ClpCopyOfArray(rhs.fromIndex_,getNumRows()+1-numberStaticRows_);
351	lowerSet_ = ClpCopyOfArray(rhs.lowerSet_,numberSets_);
352	upperSet_ = ClpCopyOfArray(rhs.upperSet_,numberSets_);
353	status_ = ClpCopyOfArray(rhs.status_,numberSets_);
354	model_ = rhs.model_;
355	sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_;
356	sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_;
357	sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_;
358	sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_;
359	numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_;
360	numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_;
361	savedBestGubDual_=rhs.savedBestGubDual_;
362	savedBestSet_=rhs.savedBestSet_;
363	noCheck_ = rhs.noCheck_;
364	infeasibilityWeight_=rhs.infeasibilityWeight_;
365	// Now secondary data
366	numberGubColumns_ = rhs.numberGubColumns_;
367	maximumGubColumns_=rhs.maximumGubColumns_;
368	maximumElements_ = rhs.maximumElements_;
369	startSet_ = ClpCopyOfArray(rhs.startSet_,numberSets_);
370	next_ = ClpCopyOfArray(rhs.next_,maximumGubColumns_);
371	startColumn_ = ClpCopyOfArray(rhs.startColumn_,maximumGubColumns_+1);
372	row_ = ClpCopyOfArray(rhs.row_,maximumElements_);
373	element_ = ClpCopyOfArray(rhs.element_,maximumElements_);
374	cost_ = ClpCopyOfArray(rhs.cost_,maximumGubColumns_);
375	id_ = ClpCopyOfArray(rhs.id_,lastDynamic_-firstDynamic_);
376	columnLower_ = ClpCopyOfArray(rhs.columnLower_,maximumGubColumns_);
377	columnUpper_ = ClpCopyOfArray(rhs.columnUpper_,maximumGubColumns_);
378	dynamicStatus_ = ClpCopyOfArray(rhs.dynamicStatus_,maximumGubColumns_);
379	}
380	return *this;
381	}
382	//-------------------------------------------------------------------
383	// Clone
384	//-------------------------------------------------------------------
385	ClpMatrixBase * ClpDynamicMatrix::clone() const
386	{
387	return new ClpDynamicMatrix(*this);
388	}
389	// Partial pricing
390	void
391	ClpDynamicMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
392	int & bestSequence, int & numberWanted)
393	{
394	numberWanted=currentWanted_;
395	assert(!model->rowScale());
396	if (numberSets_) {
397	// Do packed part before gub
398	// always???
399	//printf("normal packed price - start %d end %d (passed end %d, first %d)\n",
400	ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
401	} else {
402	// no gub
403	ClpPackedMatrix::partialPricing(model,startFraction,endFraction,bestSequence,numberWanted);
404	return;
405	}
406	if (numberWanted>0) {
407	// and do some proportion of full set
408	int startG2 = (int) (startFraction*numberSets_);
409	int endG2 = (int) (endFraction*numberSets_+0.1);
410	endG2 = CoinMin(endG2,numberSets_);
411	//printf("gub price - set start %d end %d\n",
412	// startG2,endG2);
413	double tolerance=model->currentDualTolerance();
414	double * reducedCost = model->djRegion();
415	const double * duals = model->dualRowSolution();
416	double bestDj;
417	int numberRows = model->numberRows();
418	int slackOffset = lastDynamic_+numberRows;
419	int structuralOffset = slackOffset+numberSets_;
420	// If nothing found yet can go all the way to end
421	int endAll = endG2;
422	if (bestSequence<0&&!startG2)
423	endAll = numberSets_;
424	if (bestSequence>=0) {
425	if (bestSequence!=savedBestSequence_)
426	bestDj = fabs(reducedCost[bestSequence]); // dj from slacks or permanent
427	else
428	bestDj = savedBestDj_;
429	} else {
430	bestDj=tolerance;
431	}
432	int saveSequence = bestSequence;
433	double djMod=0.0;
434	double bestDjMod=0.0;
435	//printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(),
436	// startG2,endG2,numberWanted);
437	int bestSet=-1;
438	#if 0
439	// make sure first available is clean (in case last iteration rejected)
440	cost[firstAvailable_]=0.0;
441	length[firstAvailable_]=0;
442	model->nonLinearCost()->setOne(firstAvailable_,0.0,0.0,COIN_DBL_MAX,0.0);
443	model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
444	{
445	for (int i=firstAvailable_;i<lastDynamic_;i++)
446	assert(!cost[i]);
447	}
448	#endif
449	int minSet = minimumObjectsScan_<0 ? 5 : minimumObjectsScan_;
450	int minNeg = minimumGoodReducedCosts_<0 ? 5 : minimumGoodReducedCosts_;
451	for (int iSet = startG2;iSet<endAll;iSet++) {
452	if (numberWanted+minNeg<originalWanted_&&iSet>startG2+minSet) {
453	// give up
454	numberWanted=0;
455	break;
456	} else if (iSet==endG2&&bestSequence>=0) {
457	break;
458	}
459	int gubRow = toIndex_[iSet];
460	if (gubRow>=0) {
461	djMod = duals[gubRow+numberStaticRows_]; // have I got sign right?
462	} else {
463	int iBasic = keyVariable_[iSet];
464	if (iBasic>=maximumGubColumns_) {
465	djMod = 0.0; // set not in
466	} else {
467	// get dj without
468	djMod=0.0;
469	for (CoinBigIndex j=startColumn_[iBasic];
470	j<startColumn_[iBasic+1];j++) {
471	int jRow=row_[j];
472	djMod -= duals[jRow]*element_[j];
473	}
474	djMod += cost_[iBasic];
475	// See if gub slack possible - dj is djMod
476	if (getStatus(iSet)==ClpSimplex::atLowerBound) {
477	double value = -djMod;
478	if (value>tolerance) {
479	numberWanted--;
480	if (value>bestDj) {
481	// check flagged variable and correct dj
482	if (!flagged(iSet)) {
483	bestDj=value;
484	bestSequence = slackOffset+iSet;
485	bestDjMod = djMod;
486	bestSet = iSet;
487	} else {
488	// just to make sure we don't exit before got something
489	numberWanted++;
490	abort();
491	}
492	}
493	}
494	} else if (getStatus(iSet)==ClpSimplex::atUpperBound) {
495	double value = djMod;
496	if (value>tolerance) {
497	numberWanted--;
498	if (value>bestDj) {
499	// check flagged variable and correct dj
500	if (!flagged(iSet)) {
501	bestDj=value;
502	bestSequence = slackOffset+iSet;
503	bestDjMod = djMod;
504	bestSet = iSet;
505	} else {
506	// just to make sure we don't exit before got something
507	numberWanted++;
508	abort();
509	}
510	}
511	}
512	}
513	}
514	}
515	int iSequence= startSet_[iSet];
516	while (iSequence>=0) {
517	DynamicStatus status = getDynamicStatus(iSequence);
518	if (status==atLowerBound\|\|status==atUpperBound) {
519	double value=cost_[iSequence]-djMod;
520	for (CoinBigIndex j=startColumn_[iSequence];
521	j<startColumn_[iSequence+1];j++) {
522	int jRow=row_[j];
523	value -= duals[jRow]*element_[j];
524	}
525	// change sign if at lower bound
526	if (status==atLowerBound)
527	value = -value;
528	if (value>tolerance) {
529	numberWanted--;
530	if (value>bestDj) {
531	// check flagged variable and correct dj
532	if (!flagged(iSequence)) {
533	bestDj=value;
534	bestSequence = structuralOffset+iSequence;
535	bestDjMod = djMod;
536	bestSet = iSet;
537	} else {
538	// just to make sure we don't exit before got something
539	numberWanted++;
540	}
541	}
542	}
543	}
544	iSequence = next_[iSequence]; //onto next in set
545	}
546	if (numberWanted<=0) {
547	numberWanted=0;
548	break;
549	}
550	}
551	if (bestSequence!=saveSequence) {
552	savedBestGubDual_ = bestDjMod;
553	savedBestDj_ = bestDj;
554	savedBestSequence_ = bestSequence;
555	savedBestSet_ = bestSet;
556	}
557	// See if may be finished
558	if (!startG2&&bestSequence<0)
559	infeasibilityWeight_=model_->infeasibilityCost();
560	else if (bestSequence>=0)
561	infeasibilityWeight_=-1.0;
562	}
563	currentWanted_=numberWanted;
564	}
565	/* Returns effective RHS if it is being used. This is used for long problems
566	or big gub or anywhere where going through full columns is
567	expensive. This may re-compute */
568	double *
569	ClpDynamicMatrix::rhsOffset(ClpSimplex * model,bool forceRefresh,
570	bool check)
571	{
572	//forceRefresh=true;
573	if (!model_->numberIterations())
574	forceRefresh=true;
575	//check=false;
576	#ifdef CLP_DEBUG
577	double * saveE=NULL;
578	if (rhsOffset_&&check) {
579	int numberRows = model->numberRows();
580	saveE = new double[numberRows];
581	}
582	#endif
583	if (rhsOffset_) {
584	#ifdef CLP_DEBUG
585	if (check) {
586	// no need - but check anyway
587	int numberRows = model->numberRows();
588	double * rhs = new double[numberRows];
589	int iRow;
590	int iSet;
591	CoinZeroN(rhs,numberRows);
592	// do ones at bounds before gub
593	const double * smallSolution = model->solutionRegion();
594	const double * element =matrix_->getElements();
595	const int * row = matrix_->getIndices();
596	const CoinBigIndex * startColumn = matrix_->getVectorStarts();
597	const int * length = matrix_->getVectorLengths();
598	int iColumn;
599	double objectiveOffset=0.0;
600	for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
601	if (model->getStatus(iColumn)!=ClpSimplex::basic) {
602	double value = smallSolution[iColumn];
603	for (CoinBigIndex j=startColumn[iColumn];
604	j<startColumn[iColumn]+length[iColumn];j++) {
605	int jRow=row[j];
606	rhs[jRow] -= value*element[j];
607	}
608	}
609	}
610	if (columnLower_\|\|columnUpper_) {
611	double * solution = new double [numberGubColumns_];
612	for (iSet=0;iSet<numberSets_;iSet++) {
613	int j= startSet_[iSet];
614	while (j>=0) {
615	double value=0.0;
616	if (getDynamicStatus(j)!=inSmall) {
617	if (getDynamicStatus(j)==atLowerBound) {
618	if (columnLower_)
619	value= columnLower_[j];
620	} else if (getDynamicStatus(j)==atUpperBound) {
621	value= columnUpper_[j];
622	} else if (getDynamicStatus(j)==soloKey) {
623	value =keyValue(iSet);
624	}
625	objectiveOffset += value*cost_[j];
626	}
627	solution[j]=value;
628	j = next_[j]; //onto next in set
629	}
630	}
631	// ones in gub and in small problem
632	for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
633	if (model_->getStatus(iColumn)!=ClpSimplex::basic) {
634	int jFull = id_[iColumn-firstDynamic_];
635	solution[jFull]=smallSolution[iColumn];
636	}
637	}
638	for (iSet=0;iSet<numberSets_;iSet++) {
639	int kRow = toIndex_[iSet];
640	if (kRow>=0)
641	kRow += numberStaticRows_;
642	int j= startSet_[iSet];
643	while (j>=0) {
644	double value = solution[j];
645	if (value) {
646	for (CoinBigIndex k= startColumn_[j];k<startColumn_[j+1];k++) {
647	int iRow = row_[k];
648	rhs[iRow] -= element_[k]*value;
649	}
650	if (kRow>=0)
651	rhs[kRow] -= value;
652	}
653	j = next_[j]; //onto next in set
654	}
655	}
656	delete [] solution;
657	} else {
658	// bounds
659	ClpSimplex::Status iStatus;
660	for (int iSet=0;iSet<numberSets_;iSet++) {
661	int kRow = toIndex_[iSet];
662	if (kRow<0) {
663	int iColumn = keyVariable_[iSet];
664	if (iColumn<maximumGubColumns_) {
665	// key is not treated as basic
666	double b=0.0;
667	// key is structural - where is slack
668	iStatus = getStatus(iSet);
669	assert (iStatus!=ClpSimplex::basic);
670	if (iStatus==ClpSimplex::atLowerBound)
671	b=lowerSet_[iSet];
672	else
673	b=upperSet_[iSet];
674	if (b) {
675	objectiveOffset += b*cost_[iColumn];
676	for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
677	int iRow = row_[j];
678	rhs[iRow] -= element_[j]*b;
679	}
680	}
681	}
682	}
683	}
684	}
685	if (fabs(model->objectiveOffset()-objectiveOffset_-objectiveOffset)>1.0e-1)
686	printf("old offset %g, true %g\n",model->objectiveOffset(),
687	objectiveOffset_-objectiveOffset);
688	for (iRow=0;iRow<numberRows;iRow++) {
689	if (fabs(rhs[iRow]-rhsOffset_[iRow])>1.0e-3)
690	printf("** bad effective %d - true %g old %g\n",iRow,rhs[iRow],rhsOffset_[iRow]);
691	}
692	memcpy(saveE,rhs,numberRows*sizeof(double));
693	delete [] rhs;
694	}
695	#endif
696	if (forceRefresh\|\|(refreshFrequency_&&model->numberIterations()>=
697	lastRefresh_+refreshFrequency_)) {
698	int numberRows = model->numberRows();
699	int iSet;
700	CoinZeroN(rhsOffset_,numberRows);
701	// do ones at bounds before gub
702	const double * smallSolution = model->solutionRegion();
703	const double * element =matrix_->getElements();
704	const int * row = matrix_->getIndices();
705	const CoinBigIndex * startColumn = matrix_->getVectorStarts();
706	const int * length = matrix_->getVectorLengths();
707	int iColumn;
708	double objectiveOffset=0.0;
709	for (iColumn=0;iColumn<firstDynamic_;iColumn++) {
710	if (model->getStatus(iColumn)!=ClpSimplex::basic) {
711	double value = smallSolution[iColumn];
712	for (CoinBigIndex j=startColumn[iColumn];
713	j<startColumn[iColumn]+length[iColumn];j++) {
714	int jRow=row[j];
715	rhsOffset_[jRow] -= value*element[j];
716	}
717	}
718	}
719	if (columnLower_\|\|columnUpper_) {
720	double * solution = new double [numberGubColumns_];
721	for (iSet=0;iSet<numberSets_;iSet++) {
722	int j= startSet_[iSet];
723	while (j>=0) {
724	double value=0.0;
725	if (getDynamicStatus(j)!=inSmall) {
726	if (getDynamicStatus(j)==atLowerBound) {
727	if (columnLower_)
728	value= columnLower_[j];
729	} else if (getDynamicStatus(j)==atUpperBound) {
730	value= columnUpper_[j];
731	} else if (getDynamicStatus(j)==soloKey) {
732	value =keyValue(iSet);
733	}
734	objectiveOffset += value*cost_[j];
735	}
736	solution[j]=value;
737	j = next_[j]; //onto next in set
738	}
739	}
740	// ones in gub and in small problem
741	for (iColumn=firstDynamic_;iColumn<firstAvailable_;iColumn++) {
742	if (model_->getStatus(iColumn)!=ClpSimplex::basic) {
743	int jFull = id_[iColumn-firstDynamic_];
744	solution[jFull]=smallSolution[iColumn];
745	}
746	}
747	for (iSet=0;iSet<numberSets_;iSet++) {
748	int kRow = toIndex_[iSet];
749	if (kRow>=0)
750	kRow += numberStaticRows_;
751	int j= startSet_[iSet];
752	while (j>=0) {
753	double value = solution[j];
754	if (value) {
755	for (CoinBigIndex k= startColumn_[j];k<startColumn_[j+1];k++) {
756	int iRow = row_[k];
757	rhsOffset_[iRow] -= element_[k]*value;
758	}
759	if (kRow>=0)
760	rhsOffset_[kRow] -= value;
761	}
762	j = next_[j]; //onto next in set
763	}
764	}
765	delete [] solution;
766	} else {
767	// bounds
768	ClpSimplex::Status iStatus;
769	for (int iSet=0;iSet<numberSets_;iSet++) {
770	int kRow = toIndex_[iSet];
771	if (kRow<0) {
772	int iColumn = keyVariable_[iSet];
773	if (iColumn<maximumGubColumns_) {
774	// key is not treated as basic
775	double b=0.0;
776	// key is structural - where is slack
777	iStatus = getStatus(iSet);
778	assert (iStatus!=ClpSimplex::basic);
779	if (iStatus==ClpSimplex::atLowerBound)
780	b=lowerSet_[iSet];
781	else
782	b=upperSet_[iSet];
783	if (b) {
784	objectiveOffset += b*cost_[iColumn];
785	for (CoinBigIndex j= startColumn_[iColumn];j<startColumn_[iColumn+1];j++) {
786	int iRow = row_[j];
787	rhsOffset_[iRow] -= element_[j]*b;
788	}
789	}
790	}
791	}
792	}
793	}
794	model->setObjectiveOffset(objectiveOffset_-objectiveOffset);
795	#ifdef CLP_DEBUG
796	if (saveE) {
797	int iRow;
798	for (iRow=0;iRow<numberRows;iRow++) {
799	if (fabs(saveE[iRow]-rhsOffset_[iRow])>1.0e-3)
800	printf("** %d - old eff %g new %g\n",iRow,saveE[iRow],rhsOffset_[iRow]);
801	}
802	delete [] saveE;
803	}
804	#endif
805	lastRefresh_ = model->numberIterations();
806	}
807	}
808	return rhsOffset_;
809	}
810	/*
811	update information for a pivot (and effective rhs)
812	*/
813	int
814	ClpDynamicMatrix::updatePivot(ClpSimplex * model,double oldInValue, double oldOutValue)
815	{
816
817	// now update working model
818	int sequenceIn = model->sequenceIn();
819	int sequenceOut = model->sequenceOut();
820	int numberColumns = model->numberColumns();
821	if (sequenceIn!=sequenceOut&&sequenceIn<numberColumns)
822	backToPivotRow_[sequenceIn]=model->pivotRow();
823	if (sequenceIn>=firstDynamic_&&sequenceIn<numberColumns) {
824	int bigSequence = id_[sequenceIn-firstDynamic_];
825	if (getDynamicStatus(bigSequence)!=inSmall) {
826	firstAvailable_++;
827	setDynamicStatus(bigSequence,inSmall);
828	}
829	}
830	// make sure slack is synchronized
831	if (sequenceIn>=numberColumns+numberStaticRows_) {
832	int iDynamic = sequenceIn-numberColumns-numberStaticRows_;
833	int iSet = fromIndex_[iDynamic];
834	setStatus(iSet,model->getStatus(sequenceIn));
835	}
836	if (sequenceOut>=numberColumns+numberStaticRows_) {
837	int iDynamic = sequenceOut-numberColumns-numberStaticRows_;
838	int iSet = fromIndex_[iDynamic];
839	// out may have gone through barrier - so check
840	double valueOut = model->lowerRegion()[sequenceOut];
841	if (fabs(valueOut -(double) lowerSet_[iSet])<
842	fabs(valueOut -(double) upperSet_[iSet]))
843	setStatus(iSet,ClpSimplex::atLowerBound);
844	else
845	setStatus(iSet,ClpSimplex::atUpperBound);
846	if (lowerSet_[iSet]==upperSet_[iSet])
847	setStatus(iSet,ClpSimplex::isFixed);
848	if (getStatus(iSet)!=model->getStatus(sequenceOut))
849	printf("** set %d status %d, var status %d\n",iSet,
850	getStatus(iSet),model->getStatus(sequenceOut));
851	}
852	ClpMatrixBase::updatePivot(model,oldInValue,oldOutValue);
853	#ifdef CLP_DEBUG
854	char * inSmall = new char [numberGubColumns_];
855	memset(inSmall,0,numberGubColumns_);
856	const double * solution = model->solutionRegion();
857	for (int i=0;i<numberGubColumns_;i++)
858	if (getDynamicStatus(i)==ClpDynamicMatrix::inSmall)
859	inSmall[i]=1;
860	for (int i=firstDynamic_;i<firstAvailable_;i++) {
861	int k=id_[i-firstDynamic_];
862	inSmall[k]=0;
863	//if (k>=23289&&k<23357&&solution[i])
864	//printf("var %d (in small %d) has value %g\n",k,i,solution[i]);
865	}
866	for (int i=0;i<numberGubColumns_;i++)
867	assert (!inSmall[i]);
868	delete [] inSmall;
869	#ifndef NDEBUG
870	for (int i=0;i<numberActiveSets_;i++) {
871	int iSet = fromIndex_[i];
872	assert (toIndex_[iSet]==i);
873	//if (getStatus(iSet)!=model->getRowStatus(i+numberStaticRows_))
874	//printf("*** set %d status %d, var status %d\n",iSet,
875	// getStatus(iSet),model->getRowStatus(i+numberStaticRows_));
876	//assert (model->getRowStatus(i+numberStaticRows_)==getStatus(iSet));
877	//if (iSet==1035) {
878	//printf("rhs for set %d (%d) is %g %g - cost %g\n",iSet,i,model->lowerRegion(0)[i+numberStaticRows_],
879	// model->upperRegion(0)[i+numberStaticRows_],model->costRegion(0)[i+numberStaticRows_]);
880	//}
881	}
882	#endif
883	#endif
884	return 0;
885	}
886	/*
887	utility dual function for dealing with dynamic constraints
888	mode=n see ClpGubMatrix.hpp for definition
889	Remember to update here when settled down
890	*/
891	void
892	ClpDynamicMatrix::dualExpanded(ClpSimplex * model,
893	CoinIndexedVector * array,
894	double * other,int mode)
895	{
896	switch (mode) {
897	// modify costs before transposeUpdate
898	case 0:
899	break;
900	// create duals for key variables (without check on dual infeasible)
901	case 1:
902	break;
903	// as 1 but check slacks and compute djs
904	case 2:
905	{
906	// do pivots
907	int * pivotVariable = model->pivotVariable();
908	int numberRows=numberStaticRows_+numberActiveSets_;
909	int numberColumns = model->numberColumns();
910	for (int iRow=0;iRow<numberRows;iRow++) {
911	int iPivot=pivotVariable[iRow];
912	if (iPivot<numberColumns)
913	backToPivotRow_[iPivot]=iRow;
914	}
915	if (noCheck_>=0) {
916	if (infeasibilityWeight_!=model_->infeasibilityCost()) {
917	// don't bother checking
918	sumDualInfeasibilities_=100.0;
919	numberDualInfeasibilities_=1;
920	sumOfRelaxedDualInfeasibilities_ =100.0;
921	return;
922	}
923	}
924	// In theory we should subtract out ones we have done but ....
925	// If key slack then dual 0.0
926	// If not then slack could be dual infeasible
927	// dj for key is zero so that defines dual on set
928	int i;
929	double * dual = model->dualRowSolution();
930	double dualTolerance = model->dualTolerance();
931	double relaxedTolerance=dualTolerance;
932	// we can't really trust infeasibilities if there is dual error
933	double error = CoinMin(1.0e-3,model->largestDualError());
934	// allow tolerance at least slightly bigger than standard
935	relaxedTolerance = relaxedTolerance + error;
936	// but we will be using difference
937	relaxedTolerance -= dualTolerance;
938	sumDualInfeasibilities_=0.0;
939	numberDualInfeasibilities_=0;
940	sumOfRelaxedDualInfeasibilities_ =0.0;
941	for (i=0;i<numberSets_;i++) {
942	double value=0.0;
943	int gubRow = toIndex_[i];
944	if (gubRow<0) {
945	int kColumn = keyVariable_[i];
946	if (kColumn<maximumGubColumns_) {
947	// dj without set
948	value = cost_[kColumn];
949	for (CoinBigIndex j=startColumn_[kColumn];
950	j<startColumn_[kColumn+1];j++) {
951	int iRow = row_[j];
952	value -= dual[iRow]*element_[j];
953	}
954	double infeasibility = 0.0;
955	if (getStatus(i)==ClpSimplex::atLowerBound) {
956	if (-value>dualTolerance)
957	infeasibility=-value-dualTolerance;
958	} else if (getStatus(i)==ClpSimplex::atUpperBound) {
959	if (value>dualTolerance)
960	infeasibility=value-dualTolerance;
961	}
962	if (infeasibility>0.0) {
963	sumDualInfeasibilities_ += infeasibility;
964	if (infeasibility>relaxedTolerance)
965	sumOfRelaxedDualInfeasibilities_ += infeasibility;
966	numberDualInfeasibilities_ ++;
967	}
968	}
969	} else {
970	value = dual[gubRow+numberStaticRows_];
971	}
972	// Now subtract out from all
973	int k= startSet_[i];
974	while (k>=0) {
975	if (getDynamicStatus(k)!=inSmall) {
976	double djValue = cost_[k]-value;
977	for (CoinBigIndex j=startColumn_[k];
978	j<startColumn_[k+1];j++) {
979	int iRow = row_[j];
980	djValue -= dual[iRow]*element_[j];
981	}
982	double infeasibility=0.0;
983	if (getDynamicStatus(k)==atLowerBound) {
984	if (djValue<-dualTolerance)
985	infeasibility=-djValue-dualTolerance;
986	} else if (getDynamicStatus(k)==atUpperBound) {
987	// at upper bound
988	if (djValue>dualTolerance)
989	infeasibility=djValue-dualTolerance;
990	}
991	if (infeasibility>0.0) {
992	sumDualInfeasibilities_ += infeasibility;
993	if (infeasibility>relaxedTolerance)
994	sumOfRelaxedDualInfeasibilities_ += infeasibility;
995	numberDualInfeasibilities_ ++;
996	}
997	}
998	k = next_[k]; //onto next in set
999	}
1000	}
1001	}
1002	infeasibilityWeight_=-1.0;
1003	break;
1004	// Report on infeasibilities of key variables
1005	case 3:
1006	{
1007	model->setSumDualInfeasibilities(model->sumDualInfeasibilities()+
1008	sumDualInfeasibilities_);
1009	model->setNumberDualInfeasibilities(model->numberDualInfeasibilities()+
1010	numberDualInfeasibilities_);
1011	model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities()+
1012	sumOfRelaxedDualInfeasibilities_);
1013	}
1014	break;
1015	// modify costs before transposeUpdate for partial pricing
1016	case 4:
1017	break;
1018	}
1019	}
1020	/*
1021	general utility function for dealing with dynamic constraints
1022	mode=n see ClpGubMatrix.hpp for definition
1023	Remember to update here when settled down
1024	*/
1025	int
1026	ClpDynamicMatrix::generalExpanded(ClpSimplex * model,int mode,int &number)
1027	{
1028	int returnCode=0;
1029	switch (mode) {
1030	// Fill in pivotVariable
1031	case 0:
1032	{
1033	// If no effective rhs - form it
1034	if (!rhsOffset_) {
1035	rhsOffset_ = new double[model->numberRows()];
1036	rhsOffset(model,true);
1037	}
1038	int i;
1039	int numberBasic=number;
1040	int numberColumns = model->numberColumns();
1041	// Use different array so can build from true pivotVariable_
1042	//int * pivotVariable = model->pivotVariable();
1043	int * pivotVariable = model->rowArray(0)->getIndices();
1044	for (i=0;i<numberColumns;i++) {
1045	if (model->getColumnStatus(i) == ClpSimplex::basic)
1046	pivotVariable[numberBasic++]=i;
1047	}
1048	number = numberBasic;
1049	}
1050	break;
1051	// Do initial extra rows + maximum basic
1052	case 2:
1053	{
1054	number = model->numberRows();
1055	}
1056	break;
1057	// Before normal replaceColumn
1058	case 3:
1059	{
1060	if (numberActiveSets_+numberStaticRows_==model_->numberRows()) {
1061	// no space - re-factorize
1062	returnCode=4;
1063	number=-1; // say no need for normal replaceColumn
1064	}
1065	}
1066	break;
1067	// To see if can dual or primal
1068	case 4:
1069	{
1070	returnCode= 1;
1071	}
1072	break;
1073	// save status
1074	case 5:
1075	{
1076	}
1077	break;
1078	// restore status
1079	case 6:
1080	{
1081	// maybe mismatch - redo all in small model
1082	//for (int i=firstDynamic_;i<firstAvailable_;i++) {
1083	//int jColumn = id_[i-firstDynamic_];
1084	//setDynamicStatus(jColumn,inSmall);
1085	//}
1086	}
1087	break;
1088	// unflag all variables
1089	case 8:
1090	{
1091	for (int i=0;i<numberGubColumns_;i++) {
1092	if (flagged(i)) {
1093	unsetFlagged(i);
1094	returnCode++;
1095	}
1096	}
1097	}
1098	break;
1099	// redo costs in primal
1100	case 9:
1101	{
1102	double * cost = model->costRegion();
1103	double * solution = model->solutionRegion();
1104	double * columnLower = model->lowerRegion();
1105	double * columnUpper = model->upperRegion();
1106	int i;
1107	bool doCosts = (number&4)!=0;
1108	bool doBounds = (number&1)!=0;
1109	for ( i=firstDynamic_;i<firstAvailable_;i++) {
1110	int jColumn = id_[i-firstDynamic_];
1111	if (doBounds) {
1112	if (!columnLower_&&!columnUpper_) {
1113	columnLower[i] = 0.0;
1114	columnUpper[i] = COIN_DBL_MAX;
1115	} else {
1116	if (columnLower_)
1117	columnLower[i] = columnLower_[jColumn];
1118	else
1119	columnLower[i] = 0.0;
1120	if (columnUpper_)
1121	columnUpper[i] = columnUpper_[jColumn];
1122	else
1123	columnUpper[i] = COIN_DBL_MAX;
1124	}
1125	}
1126	if (doCosts) {
1127	cost[i]=cost_[jColumn];
1128	// Original bounds
1129	if (model->nonLinearCost())
1130	model->nonLinearCost()->setOne(i,solution[i],
1131	this->columnLower(jColumn),
1132	this->columnUpper(jColumn),cost_[jColumn]);
1133	}
1134	}
1135	// and active sets
1136	for (i=0;i<numberActiveSets_;i++ ) {
1137	int iSet = fromIndex_[i];
1138	int iSequence=lastDynamic_+numberStaticRows_+i;
1139	if (doBounds) {
1140	if (lowerSet_[iSet]>-1.0e20)
1141	columnLower[iSequence] = lowerSet_[iSet];
1142	else
1143	columnLower[iSequence]=-COIN_DBL_MAX;
1144	if (upperSet_[iSet]<1.0e20)
1145	columnUpper[iSequence] = upperSet_[iSet];
1146	else
1147	columnUpper[iSequence]=COIN_DBL_MAX;
1148	}
1149	if (doCosts) {
1150	if (model->nonLinearCost()) {
1151	double trueLower;
1152	if (lowerSet_[iSet]>-1.0e20)
1153	trueLower = lowerSet_[iSet];
1154	else
1155	trueLower=-COIN_DBL_MAX;
1156	double trueUpper;
1157	if (upperSet_[iSet]<1.0e20)
1158	trueUpper = upperSet_[iSet];
1159	else
1160	trueUpper=COIN_DBL_MAX;
1161	model->nonLinearCost()->setOne(iSequence,solution[iSequence],
1162	trueLower,trueUpper,0.0);
1163	}
1164	}
1165	}
1166	}
1167	break;
1168	// return 1 if there may be changing bounds on variable (column generation)
1169	case 10:
1170	{
1171	// return 1 as bounds on rhs will change
1172	returnCode = 1;
1173	}
1174	break;
1175	// make sure set is clean
1176	case 7:
1177	{
1178	// first flag
1179	if (number>=firstDynamic_&&number<lastDynamic_) {
1180	assert (number==model->sequenceIn());
1181	// id will be sitting at firstAvailable
1182	int sequence = id_[firstAvailable_-firstDynamic_];
1183	setFlagged(sequence);
1184	model->clearFlagged(firstAvailable_);
1185	}
1186	}
1187	case 11:
1188	{
1189	int sequenceIn = model->sequenceIn();
1190	if (sequenceIn>=firstDynamic_&&sequenceIn<lastDynamic_) {
1191	assert (number==model->sequenceIn());
1192	// take out variable (but leave key)
1193	double * cost = model->costRegion();
1194	double * columnLower = model->lowerRegion();
1195	double * columnUpper = model->upperRegion();
1196	double * solution = model->solutionRegion();
1197	int * length = matrix_->getMutableVectorLengths();
1198	#ifndef NDEBUG
1199	if (length[sequenceIn]) {
1200	int * row = matrix_->getMutableIndices();
1201	CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1202	int iRow = row[startColumn[sequenceIn]+length[sequenceIn]-1];
1203	int which = iRow-numberStaticRows_;
1204	assert (which>=0);
1205	int iSet = fromIndex_[which];
1206	assert (toIndex_[iSet]==which);
1207	}
1208	#endif
1209	// no need firstAvailable_--;
1210	solution[firstAvailable_]=0.0;
1211	cost[firstAvailable_]=0.0;
1212	length[firstAvailable_]=0;
1213	model->nonLinearCost()->setOne(firstAvailable_,0.0,0.0,COIN_DBL_MAX,0.0);
1214	model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
1215	columnLower[firstAvailable_]=0.0;
1216	columnUpper[firstAvailable_]=COIN_DBL_MAX;
1217
1218	// not really in small problem
1219	int iBig=id_[sequenceIn-firstDynamic_];
1220	if (model->getStatus(sequenceIn)==ClpSimplex::atLowerBound) {
1221	setDynamicStatus(iBig,atLowerBound);
1222	if (columnLower_)
1223	modifyOffset(sequenceIn,columnLower_[iBig]);
1224	} else {
1225	setDynamicStatus(iBig,atUpperBound);
1226	modifyOffset(sequenceIn,columnUpper_[iBig]);
1227	}
1228	}
1229	}
1230	break;
1231	default:
1232	break;
1233	}
1234	return returnCode;
1235	}
1236	/* Purely for column generation and similar ideas. Allows
1237	matrix and any bounds or costs to be updated (sensibly).
1238	Returns non-zero if any changes.
1239	*/
1240	int
1241	ClpDynamicMatrix::refresh(ClpSimplex * model)
1242	{
1243	// If at beginning then create initial problem
1244	if (firstDynamic_==firstAvailable_) {
1245	initialProblem();
1246	return 1;
1247	} else if (!model->nonLinearCost()) {
1248	// will be same as last time
1249	return 1;
1250	}
1251	// lookup array
1252	int * active = new int [numberActiveSets_];
1253	CoinZeroN(active,numberActiveSets_);
1254	int iColumn;
1255	int numberColumns = model->numberColumns();
1256	double * element = matrix_->getMutableElements();
1257	int * row = matrix_->getMutableIndices();
1258	CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1259	int * length = matrix_->getMutableVectorLengths();
1260	double * cost = model->costRegion();
1261	double * columnLower = model->lowerRegion();
1262	double * columnUpper = model->upperRegion();
1263	CoinBigIndex numberElements=startColumn[firstDynamic_];
1264	// first just do lookup and basic stuff
1265	int currentNumber = firstAvailable_;
1266	firstAvailable_ = firstDynamic_;
1267	double objectiveChange=0.0;
1268	double * solution = model->solutionRegion();
1269	int currentNumberActiveSets=numberActiveSets_;
1270	for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) {
1271	int iRow = row[startColumn[iColumn]+length[iColumn]-1];
1272	int which = iRow-numberStaticRows_;
1273	assert (which>=0);
1274	int iSet = fromIndex_[which];
1275	assert (toIndex_[iSet]==which);
1276	if (model->getStatus(iColumn)==ClpSimplex::basic) {
1277	active[which]++;
1278	// may as well make key
1279	keyVariable_[iSet]=id_[iColumn-firstDynamic_];
1280	}
1281	}
1282	int i;
1283	numberActiveSets_=0;
1284	int numberDeleted=0;
1285	for (i=0;i<currentNumberActiveSets;i++) {
1286	int iRow = i+numberStaticRows_;
1287	int numberActive = active[i];
1288	int iSet = fromIndex_[i];
1289	if (model->getRowStatus(iRow)==ClpSimplex::basic) {
1290	numberActive++;
1291	// may as well make key
1292	keyVariable_[iSet]=maximumGubColumns_+iSet;
1293	}
1294	if (numberActive>1) {
1295	// keep
1296	active[i]=numberActiveSets_;
1297	numberActiveSets_++;
1298	} else {
1299	active[i]=-1;
1300	}
1301	}
1302
1303	for (iColumn=firstDynamic_;iColumn<currentNumber;iColumn++) {
1304	int iRow = row[startColumn[iColumn]+length[iColumn]-1];
1305	int which = iRow-numberStaticRows_;
1306	int jColumn = id_[iColumn-firstDynamic_];
1307	if (model->getStatus(iColumn)==ClpSimplex::atLowerBound\|\|
1308	model->getStatus(iColumn)==ClpSimplex::atUpperBound) {
1309	double value = solution[iColumn];
1310	bool toLowerBound=true;
1311	if (columnUpper_) {
1312	if (!columnLower_) {
1313	if (fabs(value-columnUpper_[jColumn])<fabs(value))
1314	toLowerBound=false;
1315	} else if (fabs(value-columnUpper_[jColumn])<fabs(value-columnLower_[iColumn])) {
1316	toLowerBound=false;
1317	}
1318	}
1319	if (toLowerBound) {
1320	// throw out to lower bound
1321	if (columnLower_) {
1322	setDynamicStatus(jColumn,atLowerBound);
1323	// treat solution as if exactly at a bound
1324	double value = columnLower[iColumn];
1325	objectiveChange += cost[iColumn]*value;
1326	} else {
1327	setDynamicStatus(jColumn,atLowerBound);
1328	}
1329	} else {
1330	// throw out to upper bound
1331	setDynamicStatus(jColumn,atUpperBound);
1332	double value = columnUpper[iColumn];
1333	objectiveChange += cost[iColumn]*value;
1334	}
1335	numberDeleted++;
1336	} else {
1337	assert(model->getStatus(iColumn)!=ClpSimplex::superBasic); // deal with later
1338	int iPut = active[which];
1339	if (iPut>=0) {
1340	// move
1341	id_[firstAvailable_-firstDynamic_] = jColumn;
1342	int numberThis = startColumn_[jColumn+1]-startColumn_[jColumn]+1;
1343	length[firstAvailable_]=numberThis;
1344	cost[firstAvailable_]=cost[iColumn];
1345	columnLower[firstAvailable_]=columnLower[iColumn];
1346	columnUpper[firstAvailable_]=columnUpper[iColumn];
1347	model->nonLinearCost()->setOne(firstAvailable_,solution[iColumn],0.0,COIN_DBL_MAX,
1348	cost_[jColumn]);
1349	CoinBigIndex base = startColumn_[jColumn];
1350	for (int j=0;j<numberThis-1;j++) {
1351	row[numberElements]=row_[base+j];
1352	element[numberElements++]=element_[base+j];
1353	}
1354	row[numberElements]=iPut+numberStaticRows_;
1355	element[numberElements++]=1.0;
1356	model->setStatus(firstAvailable_,model->getStatus(iColumn));
1357	solution[firstAvailable_] = solution[iColumn];
1358	firstAvailable_++;
1359	startColumn[firstAvailable_]=numberElements;
1360	}
1361	}
1362	}
1363	// zero solution
1364	CoinZeroN(solution+firstAvailable_,currentNumber-firstAvailable_);
1365	// zero cost
1366	CoinZeroN(cost+firstAvailable_,currentNumber-firstAvailable_);
1367	// zero lengths
1368	CoinZeroN(length+firstAvailable_,currentNumber-firstAvailable_);
1369	for ( iColumn=firstAvailable_;iColumn<currentNumber;iColumn++) {
1370	model->nonLinearCost()->setOne(iColumn,0.0,0.0,COIN_DBL_MAX,0.0);
1371	model->setStatus(iColumn,ClpSimplex::atLowerBound);
1372	columnLower[iColumn]=0.0;
1373	columnUpper[iColumn]=COIN_DBL_MAX;
1374	}
1375	// move rhs and set rest to infinite
1376	numberActiveSets_=0;
1377	for (i=0;i<currentNumberActiveSets;i++) {
1378	int iSet = fromIndex_[i];
1379	assert (toIndex_[iSet]==i);
1380	int iRow = i+numberStaticRows_;
1381	int iPut=active[i];
1382	if (iPut>=0) {
1383	int in = iRow+numberColumns;
1384	int out = iPut+numberColumns+numberStaticRows_;
1385	solution[out]=solution[in];
1386	columnLower[out]=columnLower[in];
1387	columnUpper[out]=columnUpper[in];
1388	cost[out]=cost[in];
1389	double trueLower;
1390	if (lowerSet_[iSet]>-1.0e20)
1391	trueLower = lowerSet_[iSet];
1392	else
1393	trueLower=-COIN_DBL_MAX;
1394	double trueUpper;
1395	if (upperSet_[iSet]<1.0e20)
1396	trueUpper = upperSet_[iSet];
1397	else
1398	trueUpper=COIN_DBL_MAX;
1399	model->nonLinearCost()->setOne(out,solution[out],
1400	trueLower,trueUpper,0.0);
1401	model->setStatus(out,model->getStatus(in));
1402	toIndex_[iSet]=numberActiveSets_;
1403	rhsOffset_[numberActiveSets_+numberStaticRows_]= rhsOffset_[i+numberStaticRows_];
1404	fromIndex_[numberActiveSets_++]=fromIndex_[i];
1405	} else {
1406	// adjust offset
1407	// put out as key
1408	int jColumn = keyVariable_[iSet];
1409	toIndex_[iSet]=-1;
1410	if (jColumn<maximumGubColumns_) {
1411	setDynamicStatus(jColumn,soloKey);
1412	double value = keyValue(iSet);
1413	objectiveChange += cost_[jColumn]*value;
1414	modifyOffset(jColumn,-value);
1415	}
1416	}
1417	}
1418	model->setObjectiveOffset(model->objectiveOffset()-objectiveChange);
1419	delete [] active;
1420	for (i=numberActiveSets_;i<currentNumberActiveSets;i++) {
1421	int iSequence = i+numberStaticRows_+numberColumns;
1422	solution[iSequence]=0.0;
1423	columnLower[iSequence]=-COIN_DBL_MAX;
1424	columnUpper[iSequence]=COIN_DBL_MAX;
1425	cost[iSequence]=0.0;
1426	model->nonLinearCost()->setOne(iSequence,solution[iSequence],
1427	columnLower[iSequence],
1428	columnUpper[iSequence],0.0);
1429	model->setStatus(iSequence,ClpSimplex::basic);
1430	rhsOffset_[i+numberStaticRows_]=0.0;
1431	}
1432	if (currentNumberActiveSets!=numberActiveSets_\|\|numberDeleted) {
1433	// deal with pivotVariable
1434	int * pivotVariable = model->pivotVariable();
1435	int sequence=firstDynamic_;
1436	int iRow=0;
1437	int base1 = firstDynamic_;
1438	int base2 = lastDynamic_;
1439	int base3 = numberColumns+numberStaticRows_;
1440	int numberRows = numberStaticRows_+currentNumberActiveSets;
1441	while (sequence<firstAvailable_) {
1442	int iPivot=pivotVariable[iRow];
1443	while (iPivot<base1\|\|(iPivot>=base2&&iPivot<base3)) {
1444	iPivot= pivotVariable[++iRow];
1445	}
1446	pivotVariable[iRow++]=sequence;
1447	sequence++;
1448	}
1449	// move normal basic ones down
1450	int iPut=iRow;
1451	for (;iRow<numberRows;iRow++) {
1452	int iPivot=pivotVariable[iRow];
1453	if (iPivot<base1\|\|(iPivot>=base2&&iPivot<base3)) {
1454	pivotVariable[iPut++]=iPivot;
1455	}
1456	}
1457	// and basic others
1458	for (i=0;i<numberActiveSets_;i++) {
1459	if (model->getRowStatus(i+numberStaticRows_)==ClpSimplex::basic) {
1460	pivotVariable[iPut++]=i+base3;
1461	}
1462	}
1463	for (i=numberActiveSets_;i<currentNumberActiveSets;i++) {
1464	pivotVariable[iPut++]=i+base3;
1465	}
1466	assert (iPut==numberRows);
1467	}
1468	#ifdef CLP_DEBUG
1469	#if 0
1470	printf("row for set 244 is %d, row status %d value %g ",toIndex_[244],status_[244],
1471	keyValue(244));
1472	int jj=startSet_[244];
1473	while (jj>=0) {
1474	printf("var %d status %d ",jj,dynamicStatus_[jj]);
1475	jj=next_[jj];
1476	}
1477	printf("\n");
1478	#endif
1479	#ifdef CLP_DEBUG
1480	{
1481	// debug coding to analyze set
1482	int i;
1483	int kSet=-6;
1484	if (kSet>=0) {
1485	int * back = new int [numberGubColumns_];
1486	for (i=0;i<numberGubColumns_;i++)
1487	back[i]=-1;
1488	for (i=firstDynamic_;i<firstAvailable_;i++)
1489	back[id_[i-firstDynamic_]]= i;
1490	int sequence = startSet_[kSet];
1491	if (toIndex_[kSet]<0){
1492	printf("Not in - Set %d - slack status %d, key %d\n",kSet,status_[kSet],keyVariable_[kSet]);
1493	while (sequence>=0) {
1494	printf("( %d status %d ) ",sequence,getDynamicStatus(sequence));
1495	sequence = next_[sequence];
1496	}
1497	} else {
1498	int iRow=numberStaticRows_+toIndex_[kSet];
1499	printf("In - Set %d - slack status %d, key %d offset %g slack %g\n",kSet,status_[kSet],keyVariable_[kSet],
1500	rhsOffset_[iRow],model->solutionRegion(0)[iRow]);
1501	while (sequence>=0) {
1502	int iBack = back[sequence];
1503	printf("( %d status %d value %g) ",sequence,getDynamicStatus(sequence),model->solutionRegion()[iBack]);
1504	sequence = next_[sequence];
1505	}
1506	}
1507	printf("\n");
1508	delete [] back;
1509	}
1510	}
1511	#endif
1512	int n=numberActiveSets_;
1513	for (i=0;i<numberSets_;i++) {
1514	if (toIndex_[i]<0) {
1515	//assert(keyValue(i)>=lowerSet_[i]&&keyValue(i)<=upperSet_[i]);
1516	n++;
1517	}
1518	}
1519	assert (n==numberSets_);
1520	#endif
1521	return 1;
1522	}
1523	void
1524	ClpDynamicMatrix::times(double scalar,
1525	const double * x, double * y) const
1526	{
1527	if (model_->specialOptions()!=16) {
1528	ClpPackedMatrix::times(scalar,x,y);
1529	} else {
1530	int iRow;
1531	const double * element = matrix_->getElements();
1532	const int * row = matrix_->getIndices();
1533	const CoinBigIndex * startColumn = matrix_->getVectorStarts();
1534	const int * length = matrix_->getVectorLengths();
1535	int * pivotVariable = model_->pivotVariable();
1536	for (iRow=0;iRow<numberStaticRows_+numberActiveSets_;iRow++) {
1537	y[iRow] -= scalar*rhsOffset_[iRow];
1538	int iColumn = pivotVariable[iRow];
1539	if (iColumn<lastDynamic_) {
1540	CoinBigIndex j;
1541	double value = scalar*x[iColumn];
1542	if (value) {
1543	for (j=startColumn[iColumn];
1544	j<startColumn[iColumn]+length[iColumn];j++) {
1545	int jRow=row[j];
1546	y[jRow] += value*element[j];
1547	}
1548	}
1549	}
1550	}
1551	}
1552	}
1553	// Modifies rhs offset
1554	void
1555	ClpDynamicMatrix::modifyOffset(int sequence, double amount)
1556	{
1557	if (amount) {
1558	assert (rhsOffset_);
1559	CoinBigIndex j;
1560	for (j=startColumn_[sequence];j<startColumn_[sequence+1];j++) {
1561	int iRow = row_[j];
1562	rhsOffset_[iRow] += amount*element_[j];
1563	}
1564	}
1565	}
1566	// Gets key value when none in small
1567	double
1568	ClpDynamicMatrix::keyValue(int iSet) const
1569	{
1570	double value=0.0;
1571	if (toIndex_[iSet]<0) {
1572	int key = keyVariable_[iSet];
1573	if (key<maximumGubColumns_) {
1574	if (getStatus(iSet)==ClpSimplex::atLowerBound)
1575	value=lowerSet_[iSet];
1576	else
1577	value=upperSet_[iSet];
1578	int numberKey=0;
1579	int j= startSet_[iSet];
1580	while (j>=0) {
1581	DynamicStatus status = getDynamicStatus(j);
1582	assert (status!=inSmall);
1583	if (status==soloKey) {
1584	numberKey++;
1585	} else if (status==atUpperBound) {
1586	value -= columnUpper_[j];
1587	} else if (columnLower_) {
1588	value -= columnLower_[j];
1589	}
1590	j = next_[j]; //onto next in set
1591	}
1592	assert (numberKey==1);
1593	} else {
1594	int j= startSet_[iSet];
1595	while (j>=0) {
1596	DynamicStatus status = getDynamicStatus(j);
1597	assert (status!=inSmall);
1598	assert (status!=soloKey);
1599	if (status==atUpperBound) {
1600	value += columnUpper_[j];
1601	} else if (columnLower_) {
1602	value += columnLower_[j];
1603	}
1604	j = next_[j]; //onto next in set
1605	}
1606	}
1607	}
1608	return value;
1609	}
1610	// Switches off dj checking each factorization (for BIG models)
1611	void
1612	ClpDynamicMatrix::switchOffCheck()
1613	{
1614	noCheck_=0;
1615	infeasibilityWeight_=0.0;
1616	}
1617	/* Creates a variable. This is called after partial pricing and may modify matrix.
1618	May update bestSequence.
1619	*/
1620	void
1621	ClpDynamicMatrix::createVariable(ClpSimplex * model, int & bestSequence)
1622	{
1623	int numberRows = model->numberRows();
1624	int slackOffset = lastDynamic_+numberRows;
1625	int structuralOffset = slackOffset+numberSets_;
1626	int bestSequence2=savedBestSequence_-structuralOffset;
1627	if (bestSequence>=slackOffset) {
1628	double * columnLower = model->lowerRegion();
1629	double * columnUpper = model->upperRegion();
1630	double * solution = model->solutionRegion();
1631	double * reducedCost = model->djRegion();
1632	const double * duals = model->dualRowSolution();
1633	if (toIndex_[savedBestSet_]<0) {
1634	// need to put key into basis
1635	int newRow = numberActiveSets_+numberStaticRows_;
1636	model->dualRowSolution()[newRow]=savedBestGubDual_;
1637	double valueOfKey = keyValue(savedBestSet_); // done before toIndex_ set
1638	toIndex_[savedBestSet_]=numberActiveSets_;
1639	fromIndex_[numberActiveSets_++]=savedBestSet_;
1640	int iSequence=lastDynamic_+newRow;
1641	// we need to get lower and upper correct
1642	double shift=0.0;
1643	int j= startSet_[savedBestSet_];
1644	while (j>=0) {
1645	if (getDynamicStatus(j)==atUpperBound)
1646	shift += columnUpper_[j];
1647	else if (getDynamicStatus(j)==atLowerBound&&columnLower_)
1648	shift += columnLower_[j];
1649	j = next_[j]; //onto next in set
1650	}
1651	if (lowerSet_[savedBestSet_]>-1.0e20)
1652	columnLower[iSequence] = lowerSet_[savedBestSet_];
1653	else
1654	columnLower[iSequence]=-COIN_DBL_MAX;
1655	if (upperSet_[savedBestSet_]<1.0e20)
1656	columnUpper[iSequence] = upperSet_[savedBestSet_];
1657	else
1658	columnUpper[iSequence]=COIN_DBL_MAX;
1659	#ifdef CLP_DEBUG
1660	if (model_->logLevel()==63) {
1661	printf("%d in in set %d, key is %d rhs %g %g - keyvalue %g\n",
1662	bestSequence2,savedBestSet_,keyVariable_[savedBestSet_],
1663	columnLower[iSequence],columnUpper[iSequence],valueOfKey);
1664	int j= startSet_[savedBestSet_];
1665	while (j>=0) {
1666	if (getDynamicStatus(j)==atUpperBound)
1667	printf("%d atup ",j);
1668	else if (getDynamicStatus(j)==atLowerBound)
1669	printf("%d atlo ",j);
1670	else if (getDynamicStatus(j)==soloKey)
1671	printf("%d solo ",j);
1672	else
1673	abort();
1674	j = next_[j]; //onto next in set
1675	}
1676	printf("\n");
1677	}
1678	#endif
1679	if (keyVariable_[savedBestSet_]<maximumGubColumns_) {
1680	// slack not key
1681	model_->pivotVariable()[newRow]=firstAvailable_;
1682	backToPivotRow_[firstAvailable_]=newRow;
1683	model->setStatus(iSequence,getStatus(savedBestSet_));
1684	model->djRegion()[iSequence] = savedBestGubDual_;
1685	if (getStatus(savedBestSet_)==ClpSimplex::atUpperBound)
1686	solution[iSequence] = columnUpper[iSequence];
1687	else
1688	solution[iSequence] = columnLower[iSequence];
1689	// create variable and pivot in
1690	int key=keyVariable_[savedBestSet_];
1691	setDynamicStatus(key,inSmall);
1692	double * element = matrix_->getMutableElements();
1693	int * row = matrix_->getMutableIndices();
1694	CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1695	int * length = matrix_->getMutableVectorLengths();
1696	CoinBigIndex numberElements = startColumn[firstAvailable_];
1697	int numberThis = startColumn_[key+1]-startColumn_[key]+1;
1698	if (numberElements+numberThis>numberElements_) {
1699	// need to redo
1700	numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
1701	matrix_->reserve(lastDynamic_,numberElements_);
1702	element = matrix_->getMutableElements();
1703	row = matrix_->getMutableIndices();
1704	// these probably okay but be safe
1705	startColumn = matrix_->getMutableVectorStarts();
1706	length = matrix_->getMutableVectorLengths();
1707	}
1708	// already set startColumn[firstAvailable_]=numberElements;
1709	length[firstAvailable_]=numberThis;
1710	model->costRegion()[firstAvailable_]=cost_[key];
1711	CoinBigIndex base = startColumn_[key];
1712	for (int j=0;j<numberThis-1;j++) {
1713	row[numberElements]=row_[base+j];
1714	element[numberElements++]=element_[base+j];
1715	}
1716	row[numberElements]=newRow;
1717	element[numberElements++]=1.0;
1718	id_[firstAvailable_-firstDynamic_]=key;
1719	model->setObjectiveOffset(model->objectiveOffset()+cost_[key]*
1720	valueOfKey);
1721	model->solutionRegion()[firstAvailable_]= valueOfKey;
1722	model->setStatus(firstAvailable_,ClpSimplex::basic);
1723	// ***** need to adjust effective rhs
1724	if (!columnLower_&&!columnUpper_) {
1725	columnLower[firstAvailable_] = 0.0;
1726	columnUpper[firstAvailable_] = COIN_DBL_MAX;
1727	} else {
1728	if (columnLower_)
1729	columnLower[firstAvailable_] = columnLower_[key];
1730	else
1731	columnLower[firstAvailable_] = 0.0;
1732	if (columnUpper_)
1733	columnUpper[firstAvailable_] = columnUpper_[key];
1734	else
1735	columnUpper[firstAvailable_] = COIN_DBL_MAX;
1736	}
1737	model->nonLinearCost()->setOne(firstAvailable_,solution[firstAvailable_],
1738	columnLower[firstAvailable_],
1739	columnUpper[firstAvailable_],cost_[key]);
1740	startColumn[firstAvailable_+1]=numberElements;
1741	reducedCost[firstAvailable_] = 0.0;
1742	modifyOffset(key,valueOfKey);
1743	rhsOffset_[newRow] = -shift; // sign?
1744	#ifdef CLP_DEBUG
1745	model->rowArray(1)->checkClear();
1746	#endif
1747	// now pivot in
1748	unpack(model,model->rowArray(1),firstAvailable_);
1749	model->factorization()->updateColumnFT(model->rowArray(2),model->rowArray(1));
1750	double alpha = model->rowArray(1)->denseVector()[newRow];
1751	int updateStatus =
1752	model->factorization()->replaceColumn(model,
1753	model->rowArray(2),
1754	model->rowArray(1),
1755	newRow, alpha);
1756	model->rowArray(1)->clear();
1757	if (updateStatus) {
1758	if (updateStatus==3) {
1759	// out of memory
1760	// increase space if not many iterations
1761	if (model->factorization()->pivots()<
1762	0.5*model->factorization()->maximumPivots()&&
1763	model->factorization()->pivots()<400)
1764	model->factorization()->areaFactor(
1765	model->factorization()->areaFactor() * 1.1);
1766	} else {
1767	printf("Bad returncode %d from replaceColumn\n",updateStatus);
1768	}
1769	bestSequence=-1;
1770	return;
1771	}
1772	// firstAvailable_ only finally updated if good pivot (in updatePivot)
1773	// otherwise it reverts to firstAvailableBefore_
1774	firstAvailable_++;
1775	} else {
1776	// slack key
1777	model->setStatus(iSequence,ClpSimplex::basic);
1778	model->djRegion()[iSequence]=0.0;
1779	solution[iSequence]=shift;
1780	rhsOffset_[newRow] = -shift; // sign?
1781	}
1782	// correct slack
1783	model->costRegion()[iSequence] = 0.0;
1784	model->nonLinearCost()->setOne(iSequence,solution[iSequence],columnLower[iSequence],
1785	columnUpper[iSequence],0.0);
1786	}
1787	if (savedBestSequence_>=structuralOffset) {
1788	// recompute dj and create
1789	double value=cost_[bestSequence2]-savedBestGubDual_;
1790	for (CoinBigIndex jBigIndex=startColumn_[bestSequence2];
1791	jBigIndex<startColumn_[bestSequence2+1];jBigIndex++) {
1792	int jRow=row_[jBigIndex];
1793	value -= duals[jRow]*element_[jBigIndex];
1794	}
1795	int gubRow = toIndex_[savedBestSet_] + numberStaticRows_;
1796	double * element = matrix_->getMutableElements();
1797	int * row = matrix_->getMutableIndices();
1798	CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
1799	int * length = matrix_->getMutableVectorLengths();
1800	CoinBigIndex numberElements = startColumn[firstAvailable_];
1801	int numberThis = startColumn_[bestSequence2+1]-startColumn_[bestSequence2]+1;
1802	if (numberElements+numberThis>numberElements_) {
1803	// need to redo
1804	numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
1805	matrix_->reserve(lastDynamic_,numberElements_);
1806	element = matrix_->getMutableElements();
1807	row = matrix_->getMutableIndices();
1808	// these probably okay but be safe
1809	startColumn = matrix_->getMutableVectorStarts();
1810	length = matrix_->getMutableVectorLengths();
1811	}
1812	// already set startColumn[firstAvailable_]=numberElements;
1813	length[firstAvailable_]=numberThis;
1814	model->costRegion()[firstAvailable_]=cost_[bestSequence2];
1815	CoinBigIndex base = startColumn_[bestSequence2];
1816	for (int j=0;j<numberThis-1;j++) {
1817	row[numberElements]=row_[base+j];
1818	element[numberElements++]=element_[base+j];
1819	}
1820	row[numberElements]=gubRow;
1821	element[numberElements++]=1.0;
1822	id_[firstAvailable_-firstDynamic_]=bestSequence2;
1823	//printf("best %d\n",bestSequence2);
1824	model->solutionRegion()[firstAvailable_]=0.0;
1825	if (!columnLower_&&!columnUpper_) {
1826	model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
1827	columnLower[firstAvailable_] = 0.0;
1828	columnUpper[firstAvailable_] = COIN_DBL_MAX;
1829	} else {
1830	DynamicStatus status = getDynamicStatus(bestSequence2);
1831	if (columnLower_)
1832	columnLower[firstAvailable_] = columnLower_[bestSequence2];
1833	else
1834	columnLower[firstAvailable_] = 0.0;
1835	if (columnUpper_)
1836	columnUpper[firstAvailable_] = columnUpper_[bestSequence2];
1837	else
1838	columnUpper[firstAvailable_] = COIN_DBL_MAX;
1839	if (status==atLowerBound) {
1840	solution[firstAvailable_]=columnLower[firstAvailable_];
1841	model->setStatus(firstAvailable_,ClpSimplex::atLowerBound);
1842	} else {
1843	solution[firstAvailable_]=columnUpper[firstAvailable_];
1844	model->setStatus(firstAvailable_,ClpSimplex::atUpperBound);
1845	}
1846	}
1847	model->setObjectiveOffset(model->objectiveOffset()+cost_[bestSequence2]*
1848	solution[firstAvailable_]);
1849	model->nonLinearCost()->setOne(firstAvailable_,solution[firstAvailable_],
1850	columnLower[firstAvailable_],
1851	columnUpper[firstAvailable_],cost_[bestSequence2]);
1852	bestSequence=firstAvailable_;
1853	// firstAvailable_ only updated if good pivot (in updatePivot)
1854	startColumn[firstAvailable_+1]=numberElements;
1855	//printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,savedBestGubDual_);
1856	reducedCost[bestSequence] = value;
1857	} else {
1858	// slack - row must just have been created
1859	assert (toIndex_[savedBestSet_]==numberActiveSets_-1);
1860	int newRow = numberStaticRows_+numberActiveSets_-1;
1861	bestSequence = lastDynamic_ + newRow;
1862	reducedCost[bestSequence] = savedBestGubDual_;
1863	}
1864	}
1865	// clear for next iteration
1866	savedBestSequence_=-1;
1867	}
1868	// Returns reduced cost of a variable
1869	double
1870	ClpDynamicMatrix::reducedCost(ClpSimplex * model,int sequence) const
1871	{
1872	int numberRows = model->numberRows();
1873	int slackOffset = lastDynamic_+numberRows;
1874	if (sequence<slackOffset)
1875	return model->djRegion()[sequence];
1876	else
1877	return savedBestDj_;
1878	}
1879	// Does gub crash
1880	void
1881	ClpDynamicMatrix::gubCrash()
1882	{
1883	// Do basis - cheapest or slack if feasible
1884	int longestSet=0;
1885	int iSet;
1886	for (iSet=0;iSet<numberSets_;iSet++) {
1887	int n=0;
1888	int j= startSet_[iSet];
1889	while (j>=0) {
1890	n++;
1891	j=next_[j];
1892	}
1893	longestSet = CoinMax(longestSet,n);
1894	}
1895	double * upper = new double[longestSet+1];
1896	double * cost = new double[longestSet+1];
1897	double * lower = new double[longestSet+1];
1898	double * solution = new double[longestSet+1];
1899	int * back = new int[longestSet+1];
1900	double tolerance = model_->primalTolerance();
1901	double objectiveOffset = 0.0;
1902	for (iSet=0;iSet<numberSets_;iSet++) {
1903	int iBasic=-1;
1904	double value=0.0;
1905	// find cheapest
1906	int numberInSet = 0;
1907	int j=startSet_[iSet];
1908	while (j>=0) {
1909	if (!columnLower_)
1910	lower[numberInSet]=0.0;
1911	else
1912	lower[numberInSet]= columnLower_[j];
1913	if (!columnUpper_)
1914	upper[numberInSet]=COIN_DBL_MAX;
1915	else
1916	upper[numberInSet]= columnUpper_[j];
1917	back[numberInSet++]=j;
1918	j = next_[j];
1919	}
1920	CoinFillN(solution,numberInSet,0.0);
1921	// and slack
1922	iBasic=numberInSet;
1923	solution[iBasic]=-value;
1924	lower[iBasic]=-upperSet_[iSet];
1925	upper[iBasic]=-lowerSet_[iSet];
1926	int kphase;
1927	if (value<lowerSet_[iSet]-tolerance\|\|value>upperSet_[iSet]+tolerance) {
1928	// infeasible
1929	kphase=0;
1930	// remember bounds are flipped so opposite to natural
1931	if (value<lowerSet_[iSet]-tolerance)
1932	cost[iBasic]=1.0;
1933	else
1934	cost[iBasic]=-1.0;
1935	CoinZeroN(cost,numberInSet);
1936	double dualTolerance =model_->dualTolerance();
1937	for (int iphase =kphase;iphase<2;iphase++) {
1938	if (iphase) {
1939	cost[numberInSet]=0.0;
1940	for (int j=0;j<numberInSet;j++)
1941	cost[j]= cost_[back[j]];
1942	}
1943	// now do one row lp
1944	bool improve=true;
1945	while (improve) {
1946	improve=false;
1947	double dual = cost[iBasic];
1948	int chosen =-1;
1949	double best=dualTolerance;
1950	int way=0;
1951	for (int i=0;i<=numberInSet;i++) {
1952	double dj = cost[i]-dual;
1953	double improvement =0.0;
1954	double distance=0.0;
1955	if (iphase\|\|i<numberInSet)
1956	assert (solution[i]>=lower[i]&&solution[i]<=upper[i]);
1957	if (dj>dualTolerance)
1958	improvement = dj*(solution[i]-lower[i]);
1959	else if (dj<-dualTolerance)
1960	improvement = dj*(solution[i]-upper[i]);
1961	if (improvement>best) {
1962	best=improvement;
1963	chosen=i;
1964	if (dj<0.0) {
1965	way = 1;
1966	distance = upper[i]-solution[i];
1967	} else {
1968	way = -1;
1969	distance = solution[i]-lower[i];
1970	}
1971	}
1972	}
1973	if (chosen>=0) {
1974	improve=true;
1975	// now see how far
1976	if (way>0) {
1977	// incoming increasing so basic decreasing
1978	// if phase 0 then go to nearest bound
1979	double distance=upper[chosen]-solution[chosen];
1980	double basicDistance;
1981	if (!iphase) {
1982	assert (iBasic==numberInSet);
1983	assert (solution[iBasic]>upper[iBasic]);
1984	basicDistance = solution[iBasic]-upper[iBasic];
1985	} else {
1986	basicDistance = solution[iBasic]-lower[iBasic];
1987	}
1988	// need extra coding for unbounded
1989	assert (CoinMin(distance,basicDistance)<1.0e20);
1990	if (distance>basicDistance) {
1991	// incoming becomes basic
1992	solution[chosen] += basicDistance;
1993	if (!iphase)
1994	solution[iBasic]=upper[iBasic];
1995	else
1996	solution[iBasic]=lower[iBasic];
1997	iBasic = chosen;
1998	} else {
1999	// flip
2000	solution[chosen]=upper[chosen];
2001	solution[iBasic] -= distance;
2002	}
2003	} else {
2004	// incoming decreasing so basic increasing
2005	// if phase 0 then go to nearest bound
2006	double distance=solution[chosen]-lower[chosen];
2007	double basicDistance;
2008	if (!iphase) {
2009	assert (iBasic==numberInSet);
2010	assert (solution[iBasic]<lower[iBasic]);
2011	basicDistance = lower[iBasic]-solution[iBasic];
2012	} else {
2013	basicDistance = upper[iBasic]-solution[iBasic];
2014	}
2015	// need extra coding for unbounded - for now just exit
2016	if (CoinMin(distance,basicDistance)>1.0e20) {
2017	printf("unbounded on set %d\n",iSet);
2018	iphase=1;
2019	iBasic=numberInSet;
2020	break;
2021	}
2022	if (distance>basicDistance) {
2023	// incoming becomes basic
2024	solution[chosen] -= basicDistance;
2025	if (!iphase)
2026	solution[iBasic]=lower[iBasic];
2027	else
2028	solution[iBasic]=upper[iBasic];
2029	iBasic = chosen;
2030	} else {
2031	// flip
2032	solution[chosen]=lower[chosen];
2033	solution[iBasic] += distance;
2034	}
2035	}
2036	if (!iphase) {
2037	if(iBasic<numberInSet)
2038	break; // feasible
2039	else if (solution[iBasic]>=lower[iBasic]&&
2040	solution[iBasic]<=upper[iBasic])
2041	break; // feasible (on flip)
2042	}
2043	}
2044	}
2045	}
2046	}
2047	// do solution i.e. bounds
2048	if (columnLower_\|\|columnUpper_) {
2049	for (int j=0;j<numberInSet;j++) {
2050	if (j!=iBasic) {
2051	objectiveOffset += solution[j]*cost[j];
2052	if (columnLower_&&columnUpper_) {
2053	if (fabs(solution[j]-columnLower_[back[j]])>
2054	fabs(solution[j]-columnUpper_[back[j]]))
2055	setDynamicStatus(back[j],atUpperBound);
2056	} else if (columnUpper_&&solution[j]>0.0) {
2057	setDynamicStatus(back[j],atUpperBound);
2058	} else {
2059	setDynamicStatus(back[j],atLowerBound);
2060	assert(!solution[j]);
2061	}
2062	}
2063	}
2064	}
2065	// convert iBasic back and do bounds
2066	if (iBasic==numberInSet) {
2067	// slack basic
2068	setStatus(iSet,ClpSimplex::basic);
2069	iBasic=iSet+maximumGubColumns_;
2070	} else {
2071	iBasic = back[iBasic];
2072	setDynamicStatus(iBasic,soloKey);
2073	// remember bounds flipped
2074	if (upper[numberInSet]==lower[numberInSet])
2075	setStatus(iSet,ClpSimplex::isFixed);
2076	else if (solution[numberInSet]==upper[numberInSet])
2077	setStatus(iSet,ClpSimplex::atLowerBound);
2078	else if (solution[numberInSet]==lower[numberInSet])
2079	setStatus(iSet,ClpSimplex::atUpperBound);
2080	else
2081	abort();
2082	}
2083	keyVariable_[iSet]=iBasic;
2084	}
2085	model_->setObjectiveOffset(objectiveOffset_-objectiveOffset);
2086	delete [] lower;
2087	delete [] solution;
2088	delete [] upper;
2089	delete [] cost;
2090	delete [] back;
2091	// make sure matrix is in good shape
2092	matrix_->orderMatrix();
2093	}
2094	// Populates initial matrix from dynamic status
2095	void
2096	ClpDynamicMatrix::initialProblem()
2097	{
2098	int iSet;
2099	double * element = matrix_->getMutableElements();
2100	int * row = matrix_->getMutableIndices();
2101	CoinBigIndex * startColumn = matrix_->getMutableVectorStarts();
2102	int * length = matrix_->getMutableVectorLengths();
2103	double * cost = model_->objective();
2104	double * solution = model_->primalColumnSolution();
2105	double * columnLower = model_->columnLower();
2106	double * columnUpper = model_->columnUpper();
2107	double * rowSolution = model_->primalRowSolution();
2108	double * rowLower = model_->rowLower();
2109	double * rowUpper = model_->rowUpper();
2110	CoinBigIndex numberElements=startColumn[firstDynamic_];
2111
2112	firstAvailable_ = firstDynamic_;
2113	numberActiveSets_=0;
2114	for (iSet=0;iSet<numberSets_;iSet++) {
2115	toIndex_[iSet]=-1;
2116	int numberActive=0;
2117	int whichKey=-1;
2118	if (getStatus(iSet)==ClpSimplex::basic)
2119	whichKey = maximumGubColumns_+iSet;
2120	else
2121	whichKey = -1;
2122	int j= startSet_[iSet];
2123	while (j>=0) {
2124	assert (getDynamicStatus(j)!=soloKey\|\|whichKey==-1);
2125	if (getDynamicStatus(j)==inSmall)
2126	numberActive++;
2127	else if (getDynamicStatus(j)==soloKey)
2128	whichKey=j;
2129	j = next_[j]; //onto next in set
2130	}
2131	if (getStatus(iSet)==ClpSimplex::basic&&numberActive)
2132	numberActive++;
2133	if (numberActive) {
2134	assert (numberActive>1);
2135	int iRow = numberActiveSets_+numberStaticRows_;
2136	rowSolution[iRow]=0.0;
2137	double lowerValue;
2138	if (lowerSet_[iSet]>-1.0e20)
2139	lowerValue = lowerSet_[iSet];
2140	else
2141	lowerValue=-COIN_DBL_MAX;
2142	double upperValue;
2143	if (upperSet_[iSet]<1.0e20)
2144	upperValue = upperSet_[iSet];
2145	else
2146	upperValue=COIN_DBL_MAX;
2147	rowLower[iRow]=lowerValue;
2148	rowUpper[iRow]=upperValue;
2149	if (getStatus(iSet)==ClpSimplex::basic) {
2150	model_->setRowStatus(iRow,ClpSimplex::basic);
2151	rowSolution[iRow]=0.0;
2152	} else if (getStatus(iSet)==ClpSimplex::atLowerBound) {
2153	model_->setRowStatus(iRow,ClpSimplex::atLowerBound);
2154	rowSolution[iRow]=lowerValue;
2155	} else {
2156	model_->setRowStatus(iRow,ClpSimplex::atUpperBound);
2157	rowSolution[iRow]=upperValue;
2158	}
2159	j= startSet_[iSet];
2160	while (j>=0) {
2161	DynamicStatus status = getDynamicStatus(j);
2162	if (status==inSmall) {
2163	int numberThis = startColumn_[j+1]-startColumn_[j]+1;
2164	if (numberElements+numberThis>numberElements_) {
2165	// need to redo
2166	numberElements_ = CoinMax(3*numberElements_/2,numberElements+numberThis);
2167	matrix_->reserve(lastDynamic_,numberElements_);
2168	element = matrix_->getMutableElements();
2169	row = matrix_->getMutableIndices();
2170	// these probably okay but be safe
2171	startColumn = matrix_->getMutableVectorStarts();
2172	length = matrix_->getMutableVectorLengths();
2173	}
2174	length[firstAvailable_]=numberThis;
2175	cost[firstAvailable_]=cost_[j];
2176	CoinBigIndex base = startColumn_[j];
2177	for (int k=0;k<numberThis-1;k++) {
2178	row[numberElements]=row_[base+k];
2179	element[numberElements++]=element_[base+k];
2180	}
2181	row[numberElements]=iRow;
2182	element[numberElements++]=1.0;
2183	id_[firstAvailable_-firstDynamic_]=j;
2184	solution[firstAvailable_]=0.0;
2185	model_->setStatus(firstAvailable_,ClpSimplex::basic);
2186	if (!columnLower_&&!columnUpper_) {
2187	columnLower[firstAvailable_] = 0.0;
2188	columnUpper[firstAvailable_] = COIN_DBL_MAX;
2189	} else {
2190	if (columnLower_)
2191	columnLower[firstAvailable_] = columnLower_[j];
2192	else
2193	columnLower[firstAvailable_] = 0.0;
2194	if (columnUpper_)
2195	columnUpper[firstAvailable_] = columnUpper_[j];
2196	else
2197	columnUpper[firstAvailable_] = COIN_DBL_MAX;
2198	if (status==atLowerBound) {
2199	solution[firstAvailable_]=columnLower[firstAvailable_];
2200	} else {
2201	solution[firstAvailable_]=columnUpper[firstAvailable_];
2202	}
2203	}
2204	firstAvailable_++;
2205	startColumn[firstAvailable_]=numberElements;
2206	}
2207	j = next_[j]; //onto next in set
2208	}
2209	model_->setRowStatus(numberActiveSets_+numberStaticRows_,getStatus(iSet));
2210	toIndex_[iSet]=numberActiveSets_;
2211	fromIndex_[numberActiveSets_++]=iSet;
2212	}
2213	assert (toIndex_[iSet]>=0\|\|whichKey>=0);
2214	keyVariable_[iSet]=whichKey;
2215	}
2216	return;
2217	}
2218	// Adds in a column to gub structure (called from descendant)
2219	int
2220	ClpDynamicMatrix::addColumn(int numberEntries,const int * row, const float * element,
2221	float cost, float lower, float upper,int iSet,
2222	DynamicStatus status)
2223	{
2224	// check if already in
2225	int j=startSet_[iSet];
2226	while (j>=0) {
2227	if (startColumn_[j+1]-startColumn_[j]==numberEntries) {
2228	const int * row2 = row_+startColumn_[j];
2229	const float * element2 = element_ + startColumn_[j];
2230	bool same=true;
2231	for (int k=0;k<numberEntries;k++) {
2232	if (row[k]!=row2[k]\|\|element[k]!=element2[k]) {
2233	same=false;
2234	break;
2235	}
2236	}
2237	if (same) {
2238	bool odd=false;
2239	if (cost!=cost_[j])
2240	odd =true;
2241	if (columnLower_&&lower!=columnLower_[j])
2242	odd=true;
2243	if (columnUpper_&&upper!=columnUpper_[j])
2244	odd=true;
2245	if (odd)
2246	printf("seems odd - same els but cost,lo,up are %g,%g,%g and %g,%g,%g\n",
2247	cost,lower,upper,cost_[j],
2248	columnLower_ ? columnLower_[j] : 0.0,
2249	columnUpper_ ? columnUpper_[j] : 1.0e100);
2250	else
2251	return j;
2252	}
2253	}
2254	j = next_[j];
2255	}
2256
2257	if (numberGubColumns_==maximumGubColumns_\|\|
2258	startColumn_[numberGubColumns_]+numberEntries>maximumElements_) {
2259	CoinBigIndex j;
2260	int i;
2261	int put=0;
2262	int numberElements=0;
2263	CoinBigIndex start=0;
2264	// compress - leave ones at ub and basic
2265	int * which = new int [numberGubColumns_];
2266	for (i=0;i<numberGubColumns_;i++) {
2267	CoinBigIndex end=startColumn_[i+1];
2268	if (getDynamicStatus(i)!=atLowerBound) {
2269	// keep in
2270	for (j=start;j<end;j++) {
2271	row_[numberElements]=row_[j];
2272	element_[numberElements++]=element_[j];
2273	}
2274	startColumn_[put+1]=numberElements;
2275	cost_[put]=cost_[i];
2276	if (columnLower_)
2277	columnLower_[put]=columnLower_[i];
2278	if (columnUpper_)
2279	columnUpper_[put]=columnUpper_[i];
2280	dynamicStatus_[put]=dynamicStatus_[i];
2281	id_[put]=id_[i];
2282	which[i]=put;
2283	put++;
2284	} else {
2285	// out
2286	which[i]=-1;
2287	}
2288	start=end;
2289	}
2290	// now redo startSet_ and next_
2291	int * newNext = new int [maximumGubColumns_];
2292	for (int jSet=0;jSet<numberSets_;jSet++) {
2293	int sequence = startSet_[jSet];
2294	while (which[sequence]<0) {
2295	// out
2296	assert (next_[sequence]>=0);
2297	sequence=next_[sequence];
2298	}
2299	startSet_[jSet]=which[sequence];
2300	int last = which[sequence];
2301	while (next_[sequence]>=0) {
2302	sequence = next_[sequence];
2303	if(which[sequence]>=0) {
2304	// keep
2305	newNext[last]=which[sequence];
2306	last=which[sequence];
2307	}
2308	}
2309	newNext[last]=-jSet-1;
2310	}
2311	delete [] next_;
2312	next_=newNext;
2313	delete [] which;
2314	abort();
2315	}
2316	CoinBigIndex start = startColumn_[numberGubColumns_];
2317	memcpy(row_+start,row,numberEntries*sizeof(int));
2318	memcpy(element_+start,element,numberEntries*sizeof(float));
2319	startColumn_[numberGubColumns_+1]=start+numberEntries;
2320	cost_[numberGubColumns_]=cost;
2321	if (columnLower_)
2322	columnLower_[numberGubColumns_]=lower;
2323	else
2324	assert (!lower);
2325	if (columnUpper_)
2326	columnUpper_[numberGubColumns_]=upper;
2327	else
2328	assert (upper>1.0e20);
2329	setDynamicStatus(numberGubColumns_,status);
2330	// Do next_
2331	j=startSet_[iSet];
2332	startSet_[iSet]=numberGubColumns_;
2333	next_[numberGubColumns_]=j;
2334	numberGubColumns_++;
2335	return numberGubColumns_-1;
2336	}
2337	// Returns which set a variable is in
2338	int
2339	ClpDynamicMatrix::whichSet (int sequence) const
2340	{
2341	while (next_[sequence]>=0)
2342	sequence = next_[sequence];
2343	int iSet = - next_[sequence]-1;
2344	return iSet;
2345	}

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

Download in other formats: