Context Navigation

source: branches/pre/ClpSimplexPrimalQuadratic.cpp @ 180

Last change on this file since 180 was 180, checked in by forrest, 18 years ago
new stuff
Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 82.4 KB

Line
1	// Copyright (C) 2003, International Business Machines
2	// Corporation and others. All Rights Reserved.
3
4	#include "CoinPragma.hpp"
5
6	#include <math.h>
7
8	#include "CoinHelperFunctions.hpp"
9	#include "ClpSimplexPrimalQuadratic.hpp"
10	#ifdef QUADRATIC
11	#include "ClpPrimalQuadraticDantzig.hpp"
12	#endif
13	#include "ClpFactorization.hpp"
14	#include "ClpNonLinearCost.hpp"
15	#include "CoinPackedMatrix.hpp"
16	#include "CoinIndexedVector.hpp"
17	#include "CoinWarmStartBasis.hpp"
18	#include "CoinMpsIO.hpp"
19	#include "ClpPrimalColumnPivot.hpp"
20	#include "ClpMessage.hpp"
21	#include <cfloat>
22	#include <cassert>
23	#include <string>
24	#include <stdio.h>
25	#include <iostream>
26	class tempMessage :
27	public CoinMessageHandler {
28
29	public:
30	virtual int print() ;
31	tempMessage(ClpSimplex * model);
32	ClpSimplex * model_;
33	};
34
35	// Constructor with pointer to model
36	tempMessage::tempMessage(ClpSimplex * model)
37	: CoinMessageHandler(),
38	model_(model)
39	{
40	}
41	int
42	tempMessage::print()
43	{
44	static int numberFeasible=0;
45	if (currentSource()=="Clp") {
46	if (currentMessage().externalNumber()==5) {
47	if (!numberFeasible&&!model_->nonLinearCost()->numberInfeasibilities()) {
48	numberFeasible++;
49	model_->setMaximumIterations(0);
50	}
51	}
52	}
53	return CoinMessageHandler::print();
54	}
55
56	// A sequential LP method
57	int
58	ClpSimplexPrimalQuadratic::primalSLP(int numberPasses, double deltaTolerance)
59	{
60	// Are we minimizing or maximizing
61	double whichWay=optimizationDirection();
62	// This is as a user would see
63
64	int numberColumns = this->numberColumns();
65	int numberRows = this->numberRows();
66	double * columnLower = this->columnLower();
67	double * columnUpper = this->columnUpper();
68	double * objective = this->objective();
69	double * solution = this->primalColumnSolution();
70
71	// Save objective
72
73	double * saveObjective = new double [numberColumns];
74	memcpy(saveObjective,objective,numberColumns*sizeof(double));
75
76	// Get list of non linear columns
77	CoinPackedMatrix * quadratic = quadraticObjective();
78	if (!quadratic) {
79	// no quadratic part
80	return primal(0);
81	}
82	int numberNonLinearColumns = 0;
83	int iColumn;
84	int * listNonLinearColumn = new int[numberColumns];
85	memset(listNonLinearColumn,0,numberColumns*sizeof(int));
86	const int * columnQuadratic = quadratic->getIndices();
87	const int * columnQuadraticStart = quadratic->getVectorStarts();
88	const int * columnQuadraticLength = quadratic->getVectorLengths();
89	const double * quadraticElement = quadratic->getElements();
90	for (iColumn=0;iColumn<numberColumns;iColumn++) {
91	int j;
92	for (j=columnQuadraticStart[iColumn];
93	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
94	int jColumn = columnQuadratic[j];
95	listNonLinearColumn[jColumn]=1;
96	listNonLinearColumn[iColumn]=1;
97	}
98	}
99	for (iColumn=0;iColumn<numberColumns;iColumn++) {
100	if(listNonLinearColumn[iColumn])
101	listNonLinearColumn[numberNonLinearColumns++]=iColumn;
102	}
103
104	if (!numberNonLinearColumns) {
105	delete [] listNonLinearColumn;
106	// no quadratic part
107	return primal(0);
108	}
109
110	// get feasible
111	if (!this->status()\|\|numberPrimalInfeasibilities())
112	primal(1);
113	// still infeasible
114	if (numberPrimalInfeasibilities())
115	return 0;
116
117	int jNon;
118	int * last[3];
119
120	double * trust = new double[numberNonLinearColumns];
121	double * trueLower = new double[numberNonLinearColumns];
122	double * trueUpper = new double[numberNonLinearColumns];
123	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
124	iColumn=listNonLinearColumn[jNon];
125	trust[jNon]=0.5;
126	trueLower[jNon]=columnLower[iColumn];
127	trueUpper[jNon]=columnUpper[iColumn];
128	if (solution[iColumn]<trueLower[jNon])
129	solution[iColumn]=trueLower[jNon];
130	else if (solution[iColumn]>trueUpper[jNon])
131	solution[iColumn]=trueUpper[jNon];
132	}
133	int iPass;
134	double lastObjective=1.0e31;
135	double * saveSolution = new double [numberColumns];
136	double * savePi = new double [numberRows];
137	unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
138	double targetDrop=1.0e31;
139	double objectiveOffset;
140	getDblParam(ClpObjOffset,objectiveOffset);
141	// 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
142	for (iPass=0;iPass<3;iPass++) {
143	last[iPass]=new int[numberNonLinearColumns];
144	for (jNon=0;jNon<numberNonLinearColumns;jNon++)
145	last[iPass][jNon]=0;
146	}
147	// goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
148	int goodMove=-2;
149	char * statusCheck = new char[numberColumns];
150	for (iPass=0;iPass<numberPasses;iPass++) {
151	// redo objective
152	double offset=0.0;
153	double objValue=-objectiveOffset;
154	double lambda=-1.0;
155	if (goodMove>=0) {
156	// get best interpolation
157	double coeff0=-objectiveOffset,coeff1=0.0,coeff2=0.0;
158	for (iColumn=0;iColumn<numberColumns;iColumn++) {
159	coeff0 += saveObjective[iColumn]*solution[iColumn];
160	coeff1 += saveObjective[iColumn]*(saveSolution[iColumn]-solution[iColumn]);
161	}
162	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
163	iColumn=listNonLinearColumn[jNon];
164	double valueI = solution[iColumn];
165	double valueISave = saveSolution[iColumn];
166	int j;
167	for (j=columnQuadraticStart[iColumn];
168	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
169	int jColumn = columnQuadratic[j];
170	double valueJ = solution[jColumn];
171	double valueJSave = saveSolution[jColumn];
172	double elementValue = 0.5*quadraticElement[j];
173	coeff0 += valueIvalueJelementValue;
174	coeff1 += (valueIvalueJSave+valueISavevalueJ-2.0valueIvalueJ)*elementValue;
175	coeff2 += (valueISavevalueJSave+valueIvalueJ-valueISavevalueJ-valueIvalueJSave)*elementValue;
176	}
177	}
178	double lambdaValue;
179	if (fabs(coeff2)<1.0e-9) {
180	if (coeff1+coeff2>=0.0)
181	lambda = 0.0;
182	else
183	lambda = 1.0;
184	} else {
185	lambda = -(0.5*coeff1)/coeff2;
186	if (lambda>1.0\|\|lambda<0.0) {
187	if (coeff1+coeff2>=0.0)
188	lambda = 0.0;
189	else
190	lambda = 1.0;
191	}
192	}
193	lambdaValue = lambdalambdacoeff2+lambda*coeff1+coeff0;
194	printf("coeffs %g %g %g - lastobj %g\n",coeff0,coeff1,coeff2,lastObjective);
195	printf("obj at saved %g, obj now %g zero deriv at %g - value %g\n",
196	coeff0+coeff1+coeff2,coeff0,lambda,lambdaValue);
197	if (!iPass) lambda=0.0;
198	if (lambda>0.0&&lambda<=1.0) {
199	// update solution
200	for (iColumn=0;iColumn<numberColumns;iColumn++)
201	solution[iColumn] = lambda * saveSolution[iColumn]
202	+ (1.0-lambda) * solution[iColumn];
203	if (lambda>0.999) {
204	memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
205	memcpy(status_,saveStatus,numberRows+numberColumns);
206	}
207	if (lambda>0.99999&&fabs(coeff1+coeff2)>1.0e-2) {
208	// tighten all
209	goodMove=-1;
210	}
211	}
212	}
213	memcpy(objective,saveObjective,numberColumns*sizeof(double));
214	for (iColumn=0;iColumn<numberColumns;iColumn++)
215	objValue += objective[iColumn]*solution[iColumn];
216	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
217	iColumn=listNonLinearColumn[jNon];
218	if (getColumnStatus(iColumn)==basic) {
219	if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
220	statusCheck[iColumn]='l';
221	else if (solution[iColumn]>columnUpper[iColumn]-1.0e-8)
222	statusCheck[iColumn]='u';
223	else
224	statusCheck[iColumn]='B';
225	} else {
226	if (solution[iColumn]<columnLower[iColumn]+1.0e-8)
227	statusCheck[iColumn]='L';
228	else
229	statusCheck[iColumn]='U';
230	}
231	double valueI = solution[iColumn];
232	int j;
233	for (j=columnQuadraticStart[iColumn];
234	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
235	int jColumn = columnQuadratic[j];
236	double valueJ = solution[jColumn];
237	double elementValue = quadraticElement[j];
238	elementValue *= 0.5;
239	objValue += valueIvalueJelementValue;
240	offset += valueIvalueJelementValue;
241	double gradientI = valueJ*elementValue;
242	double gradientJ = valueI*elementValue;
243	offset -= gradientI*valueI;
244	objective[iColumn] += gradientI;
245	offset -= gradientJ*valueJ;
246	objective[jColumn] += gradientJ;
247	}
248	}
249	printf("objective %g, objective offset %g\n",objValue,offset);
250	setDblParam(ClpObjOffset,objectiveOffset-offset);
251	objValue *= whichWay;
252	int * temp=last[2];
253	last[2]=last[1];
254	last[1]=last[0];
255	last[0]=temp;
256	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
257	iColumn=listNonLinearColumn[jNon];
258	double change = solution[iColumn]-saveSolution[iColumn];
259	if (change<-1.0e-5) {
260	if (fabs(change+trust[jNon])<1.0e-5)
261	temp[jNon]=-1;
262	else
263	temp[jNon]=-2;
264	} else if(change>1.0e-5) {
265	if (fabs(change-trust[jNon])<1.0e-5)
266	temp[jNon]=1;
267	else
268	temp[jNon]=2;
269	} else {
270	temp[jNon]=0;
271	}
272	}
273	// goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
274	double maxDelta=0.0;
275	if (goodMove>=0) {
276	if (objValue<=lastObjective)
277	goodMove=1;
278	else
279	goodMove=0;
280	} else {
281	maxDelta=1.0e10;
282	}
283	double maxGap=0.0;
284	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
285	iColumn=listNonLinearColumn[jNon];
286	maxDelta = max(maxDelta,
287	fabs(solution[iColumn]-saveSolution[iColumn]));
288	if (goodMove>0) {
289	if (last[0][jNon]*last[1][jNon]<0) {
290	// halve
291	trust[jNon] *= 0.5;
292	} else {
293	if (last[0][jNon]==last[1][jNon]&&
294	last[0][jNon]==last[2][jNon])
295	trust[jNon] *= 1.5;
296	}
297	} else if (goodMove!=-2&&trust[jNon]>10.0*deltaTolerance) {
298	trust[jNon] *= 0.2;
299	}
300	maxGap = max(maxGap,trust[jNon]);
301	}
302	std::cout<<"largest gap is "<<maxGap<<std::endl;
303	if (iPass>10000) {
304	for (jNon=0;jNon<numberNonLinearColumns;jNon++)
305	trust[jNon] *=0.0001;
306	}
307	if (goodMove>0) {
308	double drop = lastObjective-objValue;
309	std::cout<<"True drop was "<<drop<<std::endl;
310	std::cout<<"largest delta is "<<maxDelta<<std::endl;
311	if (maxDelta<deltaTolerance&&drop<1.0e-4&&goodMove&&lambda<0.99999) {
312	std::cout<<"Exiting"<<std::endl;
313	break;
314	}
315	}
316	if (!iPass)
317	goodMove=1;
318	targetDrop=0.0;
319	double * r = this->dualColumnSolution();
320	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
321	iColumn=listNonLinearColumn[jNon];
322	columnLower[iColumn]=max(solution[iColumn]
323	-trust[jNon],
324	trueLower[jNon]);
325	columnUpper[iColumn]=min(solution[iColumn]
326	+trust[jNon],
327	trueUpper[jNon]);
328	}
329	if (iPass) {
330	// get reduced costs
331	this->matrix()->transposeTimes(savePi,
332	this->dualColumnSolution());
333	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
334	iColumn=listNonLinearColumn[jNon];
335	double dj = objective[iColumn]-r[iColumn];
336	r[iColumn]=dj;
337	if (dj<0.0)
338	targetDrop -= dj*(columnUpper[iColumn]-solution[iColumn]);
339	else
340	targetDrop -= dj*(columnLower[iColumn]-solution[iColumn]);
341	}
342	} else {
343	memset(r,0,numberColumns*sizeof(double));
344	}
345	#if 0
346	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
347	iColumn=listNonLinearColumn[jNon];
348	if (statusCheck[iColumn]=='L'&&r[iColumn]<-1.0e-4) {
349	columnLower[iColumn]=max(solution[iColumn],
350	trueLower[jNon]);
351	columnUpper[iColumn]=min(solution[iColumn]
352	+trust[jNon],
353	trueUpper[jNon]);
354	} else if (statusCheck[iColumn]=='U'&&r[iColumn]>1.0e-4) {
355	columnLower[iColumn]=max(solution[iColumn]
356	-trust[jNon],
357	trueLower[jNon]);
358	columnUpper[iColumn]=min(solution[iColumn],
359	trueUpper[jNon]);
360	} else {
361	columnLower[iColumn]=max(solution[iColumn]
362	-trust[jNon],
363	trueLower[jNon]);
364	columnUpper[iColumn]=min(solution[iColumn]
365	+trust[jNon],
366	trueUpper[jNon]);
367	}
368	}
369	#endif
370	if (goodMove) {
371	memcpy(saveSolution,solution,numberColumns*sizeof(double));
372	memcpy(savePi,this->dualRowSolution(),numberRows*sizeof(double));
373	memcpy(saveStatus,status_,numberRows+numberColumns);
374
375	std::cout<<"Pass - "<<iPass
376	<<", target drop is "<<targetDrop
377	<<std::endl;
378	lastObjective = objValue;
379	if (targetDrop<1.0e-5&&goodMove&&iPass) {
380	printf("Exiting on target drop %g\n",targetDrop);
381	break;
382	}
383	{
384	double * r = this->dualColumnSolution();
385	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
386	iColumn=listNonLinearColumn[jNon];
387	printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
388	jNon,trust[jNon],iColumn,solution[iColumn],objective[iColumn],
389	r[iColumn],statusCheck[iColumn],columnLower[iColumn],
390	columnUpper[iColumn]);
391	}
392	}
393	setLogLevel(63);
394	this->scaling(false);
395	this->primal(1);
396	if (this->status()) {
397	CoinMpsIO writer;
398	writer.setMpsData(*matrix(), COIN_DBL_MAX,
399	getColLower(), getColUpper(),
400	getObjCoefficients(),
401	(const char) 0 /integrality*/,
402	getRowLower(), getRowUpper(),
403	NULL,NULL);
404	writer.writeMps("xx.mps");
405	}
406	assert (!this->status()); // must be feasible
407	goodMove=1;
408	} else {
409	// bad pass - restore solution
410	printf("Backtracking\n");
411	memcpy(solution,saveSolution,numberColumns*sizeof(double));
412	memcpy(this->dualRowSolution(),savePi,numberRows*sizeof(double));
413	memcpy(status_,saveStatus,numberRows+numberColumns);
414	iPass--;
415	goodMove=-1;
416	}
417	}
418	// restore solution
419	memcpy(solution,saveSolution,numberColumns*sizeof(double));
420	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
421	iColumn=listNonLinearColumn[jNon];
422	columnLower[iColumn]=max(solution[iColumn],
423	trueLower[jNon]);
424	columnUpper[iColumn]=min(solution[iColumn],
425	trueUpper[jNon]);
426	}
427	delete [] statusCheck;
428	delete [] savePi;
429	delete [] saveStatus;
430	// redo objective
431	double offset=0.0;
432	double objValue=-objectiveOffset;
433	memcpy(objective,saveObjective,numberColumns*sizeof(double));
434	for (iColumn=0;iColumn<numberColumns;iColumn++)
435	objValue += objective[iColumn]*solution[iColumn];
436	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
437	iColumn=listNonLinearColumn[jNon];
438	double valueI = solution[iColumn];
439	int j;
440	for (j=columnQuadraticStart[iColumn];
441	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
442	int jColumn = columnQuadratic[j];
443	double valueJ = solution[jColumn];
444	double elementValue = quadraticElement[j];
445	objValue += 0.5valueIvalueJ*elementValue;
446	offset += 0.5valueIvalueJ*elementValue;
447	double gradientI = valueJ*elementValue;
448	double gradientJ = valueI*elementValue;
449	offset -= gradientI*valueI;
450	objective[iColumn] += gradientI;
451	offset -= gradientJ*valueJ;
452	objective[jColumn] += gradientJ;
453	}
454	}
455	printf("objective %g, objective offset %g\n",objValue,offset);
456	setDblParam(ClpObjOffset,objectiveOffset-offset);
457	this->primal(1);
458	// redo values
459	setDblParam(ClpObjOffset,objectiveOffset);
460	objectiveValue_ += optimizationDirection_*offset;
461	for (jNon=0;jNon<numberNonLinearColumns;jNon++) {
462	iColumn=listNonLinearColumn[jNon];
463	columnLower[iColumn]= trueLower[jNon];
464	columnUpper[iColumn]= trueUpper[jNon];
465	}
466	memcpy(objective,saveObjective,numberColumns*sizeof(double));
467	delete [] saveSolution;
468	for (iPass=0;iPass<3;iPass++)
469	delete [] last[iPass];
470	delete [] trust;
471	delete [] trueUpper;
472	delete [] trueLower;
473	delete [] saveObjective;
474	delete [] listNonLinearColumn;
475	return 0;
476	}
477	// Dantzig's method
478	int
479	ClpSimplexPrimalQuadratic::primalQuadratic(int phase)
480	{
481	#if 0
482	// Dantzig's method looks easier - lets try that
483
484	int numberColumns = this->numberColumns();
485	double * columnLower = this->columnLower();
486	double * columnUpper = this->columnUpper();
487	double * objective = this->objective();
488	double * solution = this->primalColumnSolution();
489	double * dj = this->dualColumnSolution();
490	double * pi = this->dualRowSolution();
491
492	int numberRows = this->numberRows();
493	double * rowLower = this->rowLower();
494	double * rowUpper = this->rowUpper();
495
496	// and elements
497	CoinPackedMatrix * matrix = this->matrix();
498	const int * row = matrix->getIndices();
499	const int * columnStart = matrix->getVectorStarts();
500	const double * element = matrix->getElements();
501	const int * columnLength = matrix->getVectorLengths();
502
503	// Get list of non linear columns
504	CoinPackedMatrix * quadratic = quadraticObjective();
505	if (!quadratic\|\|!quadratic->getNumElements()) {
506	// no quadratic part
507	return primal(1);
508	}
509
510	int iColumn;
511	const int * columnQuadratic = quadratic->getIndices();
512	const int * columnQuadraticStart = quadratic->getVectorStarts();
513	const int * columnQuadraticLength = quadratic->getVectorLengths();
514	const double * quadraticElement = quadratic->getElements();
515	#if 0
516	// deliberate bad solution
517	// Change to use phase
518	//double * saveO = new double[numberColumns];
519	//memcpy(saveO,objective,numberColumns*sizeof(double));
520	//memset(objective,0,numberColumns*sizeof(double));
521	tempMessage messageHandler(this);;
522	passInMessageHandler(&messageHandler);
523	factorization()->maximumPivots(1);
524	primal();
525	CoinMessageHandler handler2;
526	passInMessageHandler(&handler2);
527	factorization()->maximumPivots(100);
528	setMaximumIterations(1000);
529	#endif
530	//memcpy(objective,saveO,numberColumns*sizeof(double));
531	//printf("For testing - deliberate bad solution\n");
532	//columnUpper[0]=0.0;
533	//columnLower[0]=0.0;
534	//quadraticSLP(50,1.0e-4);
535	//primal(1);
536	//columnUpper[0]=COIN_DBL_MAX;
537
538	// Create larger problem
539	// First workout how many rows extra
540	ClpQuadraticInfo info(this);
541	int numberQuadratic = info.numberQuadraticColumns();
542	int newNumberRows = numberRows+numberQuadratic;
543	int newNumberColumns = numberColumns + numberRows + numberQuadratic;
544	int numberElements = 2*matrix->getNumElements()
545	+2*quadratic->getNumElements()
546	+ numberQuadratic;
547	double * elements2 = new double[numberElements];
548	int * start2 = new int[newNumberColumns+1];
549	int * row2 = new int[numberElements];
550	double * objective2 = new double[newNumberColumns];
551	double * columnLower2 = new double[newNumberColumns];
552	double * columnUpper2 = new double[newNumberColumns];
553	double * rowLower2 = new double[newNumberRows];
554	double * rowUpper2 = new double[newNumberRows];
555	const int * which = info.quadraticSequence();
556	const int * back = info.backSequence();
557	memcpy(rowLower2,rowLower,numberRows*sizeof(double));
558	memcpy(rowUpper2,rowUpper,numberRows*sizeof(double));
559	int iRow;
560	for (iRow=0;iRow<numberQuadratic;iRow++) {
561	double cost = objective[back[iRow]];
562	rowLower2[iRow+numberRows]=cost;
563	rowUpper2[iRow+numberRows]=cost;
564	}
565	memset(objective2,0,newNumberColumns*sizeof(double));
566	// Get a row copy of quadratic objective in standard format
567	CoinPackedMatrix copyQ;
568	copyQ.reverseOrderedCopyOf(*quadratic);
569	const int * columnQ = copyQ.getIndices();
570	const CoinBigIndex * rowStartQ = copyQ.getVectorStarts();
571	const int * rowLengthQ = copyQ.getVectorLengths();
572	const double * elementByRowQ = copyQ.getElements();
573	// Move solution across
574	double * solution2 = new double[newNumberColumns];
575	memset(solution2,0,newNumberColumns*sizeof(double));
576	newNumberColumns=0;
577	numberElements=0;
578	start2[0]=0;
579	// x
580	memcpy(dj,objective,numberColumns*sizeof(double));
581	for (iColumn=0;iColumn<numberColumns;iColumn++) {
582	// Original matrix
583	columnLower2[iColumn]=columnLower[iColumn];
584	columnUpper2[iColumn]=columnUpper[iColumn];
585	solution2[iColumn]=solution[iColumn];
586	int j;
587	for (j=columnStart[iColumn];
588	j<columnStart[iColumn]+columnLength[iColumn];
589	j++) {
590	elements2[numberElements]=element[j];
591	row2[numberElements++]=row[j];
592	}
593	if (which[iColumn]>=0) {
594	// Quadratic and modify djs
595	for (j=columnQuadraticStart[iColumn];
596	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
597	j++) {
598	int jColumn = columnQuadratic[j];
599	double value = quadraticElement[j];
600	if (iColumn!=jColumn)
601	value *= 0.5;
602	dj[jColumn] += solution[iColumn]*value;
603	elements2[numberElements]=-value;
604	row2[numberElements++]=which[jColumn]+numberRows;
605	}
606	for (j=rowStartQ[iColumn];
607	j<rowStartQ[iColumn]+rowLengthQ[iColumn];
608	j++) {
609	int jColumn = columnQ[j];
610	double value = elementByRowQ[j];
611	if (iColumn!=jColumn) {
612	value *= 0.5;
613	dj[jColumn] += solution[iColumn]*value;
614	elements2[numberElements]=-value;
615	row2[numberElements++]=which[jColumn]+numberRows;
616	}
617	}
618	}
619	start2[iColumn+1]=numberElements;
620	}
621	newNumberColumns=numberColumns;
622	// pi
623	// Get a row copy in standard format
624	CoinPackedMatrix copy;
625	copy.reverseOrderedCopyOf(*(this->matrix()));
626	// get matrix data pointers
627	const int * column = copy.getIndices();
628	const CoinBigIndex * rowStart = copy.getVectorStarts();
629	const int * rowLength = copy.getVectorLengths();
630	const double * elementByRow = copy.getElements();
631	for (iRow=0;iRow<numberRows;iRow++) {
632	solution2[newNumberColumns]=pi[iRow];
633	double value = pi[iRow];
634	columnLower2[newNumberColumns]=-COIN_DBL_MAX;
635	columnUpper2[newNumberColumns]=COIN_DBL_MAX;
636	int j;
637	for (j=rowStart[iRow];
638	j<rowStart[iRow]+rowLength[iRow];
639	j++) {
640	double elementValue=elementByRow[j];
641	int jColumn = column[j];
642	elements2[numberElements]=elementValue;
643	row2[numberElements++]=jColumn+numberRows;
644	dj[jColumn]-= value*elementValue;
645	}
646	newNumberColumns++;
647	start2[newNumberColumns]=numberElements;
648	}
649	// djs
650	for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
651	columnLower2[newNumberColumns]=-COIN_DBL_MAX;
652	columnUpper2[newNumberColumns]=COIN_DBL_MAX;
653	solution2[newNumberColumns]=dj[iColumn];
654	elements2[numberElements]=1.0;
655	row2[numberElements++]=back[iColumn]+numberRows;
656	newNumberColumns++;
657	start2[newNumberColumns]=numberElements;
658	}
659	// Create model
660	ClpSimplex model2(*this);
661	model2.resize(0,0);
662	model2.loadProblem(newNumberColumns,newNumberRows,
663	start2,row2, elements2,
664	columnLower2,columnUpper2,
665	objective2,
666	rowLower2,rowUpper2);
667	delete [] objective2;
668	delete [] rowLower2;
669	delete [] rowUpper2;
670	delete [] columnLower2;
671	delete [] columnUpper2;
672	// Now create expanded quadratic objective for use in primalRow
673	// Later pack down in some way
674	start2[0]=0;
675	numberElements=0;
676	for (iColumn=0;iColumn<numberColumns;iColumn++) {
677	// Quadratic
678	int j;
679	for (j=columnQuadraticStart[iColumn];
680	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
681	j++) {
682	int jColumn = columnQuadratic[j];
683	double value = quadraticElement[j];
684	if (iColumn!=jColumn)
685	value *= 0.5;
686	elements2[numberElements]=value;
687	row2[numberElements++]=jColumn;
688	}
689	for (j=rowStartQ[iColumn];
690	j<rowStartQ[iColumn]+rowLengthQ[iColumn];
691	j++) {
692	int jColumn = columnQ[j];
693	double value = elementByRowQ[j];
694	if (iColumn!=jColumn) {
695	value *= 0.5;
696	elements2[numberElements]=value;
697	row2[numberElements++]=jColumn;
698	}
699	}
700	start2[iColumn+1]=numberElements;
701	}
702	// and pad
703	for (;iColumn<newNumberColumns;iColumn++)
704	start2[iColumn+1]=numberElements;
705	// Load up objective
706	model2.loadQuadraticObjective(newNumberColumns,start2,row2,elements2);
707	delete [] start2;
708	delete [] row2;
709	delete [] elements2;
710	model2.allSlackBasis();
711	model2.scaling(false);
712	model2.setLogLevel(this->logLevel());
713	// Move solution across
714	memcpy(model2.primalColumnSolution(),solution2,
715	newNumberColumns*sizeof(double));
716	columnLower2 = model2.columnLower();
717	columnUpper2 = model2.columnUpper();
718	delete [] solution2;
719	solution2 = model2.primalColumnSolution();
720	// Compute row activities and check feasible
721	double * rowSolution2 = model2.primalRowSolution();
722	memset(rowSolution2,0,newNumberRows*sizeof(double));
723	model2.times(1.0,solution2,rowSolution2);
724	rowLower2 = model2.rowLower();
725	rowUpper2 = model2.rowUpper();
726	#if 0
727	for (iColumn=0;iColumn<numberColumns;iColumn++) {
728	Status xStatus = getColumnStatus(iColumn);
729	bool isSuperBasic;
730	int iS = iColumn+newNumberRows;
731	double value = solution2[iS];
732	if (fabs(value)>dualTolerance_)
733	isSuperBasic=true;
734	else
735	isSuperBasic=false;
736	// For moment take all x out of basis
737	// Does not seem right
738	isSuperBasic=true;
739	model2.setColumnStatus(iColumn,xStatus);
740	if (xStatus==basic) {
741	if (!isSuperBasic) {
742	model2.setRowStatus(numberRows+iColumn,basic);
743	model2.setColumnStatus(iS,superBasic);
744	} else {
745	model2.setRowStatus(numberRows+iColumn,isFixed);
746	model2.setColumnStatus(iS,basic);
747	model2.setColumnStatus(iColumn,superBasic);
748	}
749	} else {
750	model2.setRowStatus(numberRows+iColumn,isFixed);
751	model2.setColumnStatus(iS,basic);
752	}
753	}
754	for (iRow=0;iRow<numberRows;iRow++) {
755	Status rowStatus = getRowStatus(iRow);
756	model2.setRowStatus(iRow,rowStatus);
757	if (rowStatus!=basic) {
758	model2.setColumnStatus(iRow+numberColumns,basic); // make dual basic
759	}
760	assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
761	assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
762	}
763	// why ?? - take duals out and adjust
764	for (iRow=0;iRow<numberRows;iRow++) {
765	model2.setRowStatus(iRow,basic);
766	model2.setColumnStatus(iRow+numberColumns,superBasic);
767	solution2[iRow+numberColumns]=0.0;
768	}
769	#else
770	for (iRow=0;iRow<numberRows;iRow++) {
771	assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
772	assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
773	model2.setRowStatus(iRow,basic);
774	model2.setColumnStatus(iRow+numberColumns,superBasic);
775	solution2[iRow+numberColumns]=0.0;
776	}
777	for (iColumn=numberRows+numberColumns;iColumn<newNumberColumns;iColumn++) {
778	model2.setColumnStatus(iColumn,basic);
779	model2.setRowStatus(iColumn-numberColumns,isFixed);
780	}
781	#endif
782	memset(rowSolution2,0,newNumberRows*sizeof(double));
783	model2.times(1.0,solution2,rowSolution2);
784	for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
785	int iS = back[iColumn]+newNumberRows;
786	int iRow = iColumn+numberRows;
787	double value = rowSolution2[iRow];
788	if (value>rowUpper2[iRow]) {
789	rowSolution2[iRow] = rowUpper2[iRow];
790	solution2[iS]-=value-rowUpper2[iRow];
791	} else {
792	rowSolution2[iRow] = rowLower2[iRow];
793	solution2[iS]-=value-rowLower2[iRow];
794	}
795	}
796
797
798	// See if any s sub j have wrong sign and/or use djs from infeasibility objective
799	double objectiveOffset;
800	getDblParam(ClpObjOffset,objectiveOffset);
801	double objValue = -objectiveOffset;
802	for (iColumn=0;iColumn<numberColumns;iColumn++)
803	objValue += objective[iColumn]*solution2[iColumn];
804	for (iColumn=0;iColumn<numberColumns;iColumn++) {
805	double valueI = solution2[iColumn];
806	int j;
807	for (j=columnQuadraticStart[iColumn];
808	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
809	int jColumn = columnQuadratic[j];
810	double valueJ = solution2[jColumn];
811	double elementValue = quadraticElement[j];
812	objValue += 0.5valueIvalueJ*elementValue;
813	}
814	}
815	printf("Objective value %g\n",objValue);
816	for (iColumn=0;iColumn<newNumberColumns;iColumn++)
817	printf("%d %g\n",iColumn,solution2[iColumn]);
818	#if 1
819	CoinMpsIO writer;
820	writer.setMpsData(*model2.matrix(), COIN_DBL_MAX,
821	model2.getColLower(), model2.getColUpper(),
822	model2.getObjCoefficients(),
823	(const char) 0 /integrality*/,
824	model2.getRowLower(), model2.getRowUpper(),
825	NULL,NULL);
826	writer.writeMps("xx.mps");
827	#endif
828	// Now do quadratic
829	// If we did not do Slp we should have primal feasible basic solution
830	// Do safe cast as no data
831	ClpSimplexPrimalQuadratic * modelPtr =
832	(ClpSimplexPrimalQuadratic *) (&model2);
833	ClpPrimalQuadraticDantzig dantzigP(modelPtr,numberRows);
834	modelPtr->setPrimalColumnPivotAlgorithm(dantzigP);
835	modelPtr->messageHandler()->setLogLevel(63);
836	modelPtr->primalQuadratic2(this,phase);
837	memcpy(dualRowSolution(),model2.primalColumnSolution()+numberColumns_,numberRows_*sizeof(double));
838	memcpy(primalColumnSolution(),model2.primalColumnSolution(),numberColumns_*sizeof(double));
839	memset(model2.primalRowSolution(),0,newNumberRows*sizeof(double));
840	model2.times(1.0,model2.primalColumnSolution(),
841	model2.primalRowSolution());
842	memcpy(dualColumnSolution(),model2.primalColumnSolution()+numberRows_+numberColumns_,
843	numberColumns_*sizeof(double));
844	objValue = -objectiveOffset;
845	for (iColumn=0;iColumn<numberColumns;iColumn++)
846	objValue += objective[iColumn]*solution2[iColumn];
847	for (iColumn=0;iColumn<numberColumns;iColumn++) {
848	double valueI = solution2[iColumn];
849	if (fabs(valueI)>1.0e-5) {
850	int djColumn = iColumn+numberRows+numberColumns;
851	assert(solution2[djColumn]<1.0e-7);
852	}
853	int j;
854	for (j=columnQuadraticStart[iColumn];
855	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
856	int jColumn = columnQuadratic[j];
857	double valueJ = solution2[jColumn];
858	double elementValue = quadraticElement[j];
859	objValue += 0.5valueIvalueJ*elementValue;
860	}
861	}
862	printf("Objective value %g\n",objValue);
863	objectiveValue_ = objValue + objectiveOffset;
864	return 0;
865	#else
866	// Get a feasible solution
867	if (numberPrimalInfeasibilities())
868	primal(1);
869	// still infeasible
870	if (numberPrimalInfeasibilities())
871	return 1;
872	ClpQuadraticInfo info;
873	ClpSimplexPrimalQuadratic * model2 = makeQuadratic(info);
874	#if 0
875	CoinMpsIO writer;
876	writer.setMpsData(*model2->matrix(), COIN_DBL_MAX,
877	model2->getColLower(), model2->getColUpper(),
878	model2->getObjCoefficients(),
879	(const char) 0 /integrality*/,
880	model2->getRowLower(), model2->getRowUpper(),
881	NULL,NULL);
882	writer.writeMps("xx.mps");
883	#endif
884	// Now do quadratic
885	ClpPrimalQuadraticDantzig dantzigP(model2,numberRows_);
886	model2->setPrimalColumnPivotAlgorithm(dantzigP);
887	model2->messageHandler()->setLogLevel(63);
888	model2->primalQuadratic2(this,phase);
889	endQuadratic(model2,info);
890	return 0;
891	#endif
892	}
893	int ClpSimplexPrimalQuadratic::primalQuadratic2 (const ClpSimplexPrimalQuadratic * originalModel,int phase)
894	{
895
896	algorithm_ = +2;
897
898	// save data
899	ClpDataSave data = saveData();
900
901	// Assume problem is feasible
902	// Stuff below will be moved into a class
903	int numberXColumns = originalModel->numberColumns();
904	int numberXRows = originalModel->numberRows();
905	int baseS = numberXColumns+numberXRows;
906	int * pair = new int [numberColumns_];
907	int iColumn;
908	for (iColumn=0;iColumn<numberXColumns;iColumn++) {
909	int iS = iColumn+baseS;
910	pair[iColumn]=iS;
911	pair[iS]=iColumn;
912	}
913	#if 0
914	// Throw out all x variables whose dj is nonzero
915	for (iColumn=0;iColumn<numberXColumns;iColumn++) {
916	if (getColumnStatus(iColumn)==basic) {
917	int jColumn = iColumn+numberRows_;
918	if (getColumnStatus(jColumn)==basic)
919	setColumnStatus(iColumn,superBasic);
920	}
921	}
922	int originalNumberRows = originalModel->numberRows();
923	for (int i=0;i<originalNumberRows;i++) {
924	int jSequence = i+numberXColumns;
925	if (getRowStatus(i)==basic)
926	setColumnStatus(jSequence,superBasic);
927	}
928	#endif
929	// Save solution
930	double * saveSolution = new double [numberRows_+numberColumns_];
931	assert (!scalingFlag_);
932	memcpy(saveSolution,columnActivity_,numberColumns_*sizeof(double));
933	memcpy(saveSolution+numberColumns_,rowActivity_,numberRows_*sizeof(double));
934	// initialize - values pass coding and algorithm_ is +1
935	if (!startup(1)) {
936
937	// Setup useful stuff
938	ClpQuadraticInfo info(originalModel);
939	if (phase)
940	memcpy(solution_,saveSolution,
941	(numberRows_+numberColumns_)*sizeof(double));
942	int lastCleaned=0; // last time objective or bounds cleaned up
943	int sequenceIn=-1;
944	int crucialSj=-1;
945
946	// special nonlinear cost
947	delete nonLinearCost_;
948	nonLinearCost_ = new ClpNonLinearCost(this,originalModel->numberColumns());
949	// Say no pivot has occurred (for steepest edge and updates)
950	pivotRow_=-2;
951
952	// This says whether to restore things etc
953	int factorType=0;
954	/*
955	Status of problem:
956	0 - optimal
957	1 - infeasible
958	2 - unbounded
959	-1 - iterating
960	-2 - factorization wanted
961	-3 - redo checking without factorization
962	-4 - looks infeasible
963	-5 - looks unbounded
964	*/
965	while (problemStatus_<0) {
966	int iRow,iColumn;
967	// clear
968	for (iRow=0;iRow<4;iRow++) {
969	rowArray_[iRow]->clear();
970	}
971
972	for (iColumn=0;iColumn<2;iColumn++) {
973	columnArray_[iColumn]->clear();
974	}
975
976	// give matrix (and model costs and bounds a chance to be
977	// refreshed (normally null)
978	matrix_->refresh(this);
979	// If we have done no iterations - special
980	if (lastGoodIteration_==numberIterations_)
981	factorType=3;
982	if (phase)
983	memcpy(saveSolution,solution_,
984	(numberRows_+numberColumns_)*sizeof(double));
985	if (phase&&0) {
986	// Clean solution
987	for (iColumn=0;iColumn<numberColumns_;iColumn++) {
988	if (getColumnStatus(iColumn)==isFree)
989	solution_[iColumn]=0.0;
990	}
991	}
992	// may factorize, checks if problem finished
993	statusOfProblemInPrimal(lastCleaned,factorType,progress_);
994	if (phase)
995	memcpy(solution_,saveSolution,
996	(numberRows_+numberColumns_)*sizeof(double));
997
998	// Compute objective function from scratch
999	const CoinPackedMatrix * quadratic = originalModel->quadraticObjective();
1000	const int * columnQuadratic = quadratic->getIndices();
1001	const int * columnQuadraticStart = quadratic->getVectorStarts();
1002	const int * columnQuadraticLength = quadratic->getVectorLengths();
1003	const double * quadraticElement = quadratic->getElements();
1004	const double * originalCost = originalModel->objective();
1005	objectiveValue_=0.0;
1006	for (iColumn=0;iColumn<numberXColumns;iColumn++) {
1007	double valueI = solution_[iColumn];
1008	objectiveValue_ += valueI*originalCost[iColumn];
1009	int j;
1010	for (j=columnQuadraticStart[iColumn];
1011	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1012	int jColumn = columnQuadratic[j];
1013	double valueJ = solution_[jColumn];
1014	double elementValue = 0.5*quadraticElement[j];
1015	objectiveValue_ += valueIvalueJelementValue;
1016	}
1017	}
1018
1019	// Say good factorization
1020	factorType=1;
1021
1022	// Say no pivot has occurred (for steepest edge and updates)
1023	pivotRow_=-2;
1024	// Check problem phase
1025	// We assume all X are feasible
1026	phase=0;
1027	if (saveSolution) {
1028	sequenceIn=-1;
1029	for (iColumn=0;iColumn<numberXColumns;iColumn++) {
1030	double value = solution_[iColumn];
1031	double lower = lower_[iColumn];
1032	double upper = upper_[iColumn];
1033	if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
1034	if (getColumnStatus(iColumn)!=basic) {
1035	if (phase!=2) {
1036	phase=2;
1037	sequenceIn=iColumn;
1038	}
1039	}
1040	}
1041	if (getColumnStatus(iColumn)==basic) {
1042	int iS = pair[iColumn];
1043	assert (iS>=0&&getColumnStatus(iS)!=basic);
1044	if(fabs(solution_[iS])>dualTolerance_) {
1045	if (phase==0) {
1046	phase=1;
1047	sequenceIn=iS;
1048	}
1049	}
1050	}
1051	}
1052	int offset=numberXColumns-numberColumns_;
1053	for (iColumn=numberColumns_;iColumn<numberColumns_+numberXRows;iColumn++) {
1054	double value = solution_[iColumn];
1055	double lower = lower_[iColumn];
1056	double upper = upper_[iColumn];
1057	if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
1058	if (getColumnStatus(iColumn)!=basic) {
1059	if (phase!=2) {
1060	phase=2;
1061	sequenceIn=iColumn;
1062	}
1063	}
1064	}
1065	if (getColumnStatus(iColumn)==basic) {
1066	int iS = iColumn+offset;
1067	assert (getColumnStatus(iS)!=basic);
1068	if(fabs(solution_[iS])>dualTolerance_) {
1069	if (phase==0) {
1070	phase=1;
1071	sequenceIn=iS;
1072	}
1073	}
1074	}
1075	}
1076	}
1077	if (!phase) {
1078	delete [] saveSolution;
1079	saveSolution=NULL;
1080	}
1081
1082	// exit if victory declared
1083	if (!phase&&primalColumnPivot_->pivotColumn(rowArray_[1],
1084	rowArray_[2],rowArray_[3],
1085	columnArray_[0],
1086	columnArray_[1]) < 0) {
1087	problemStatus_=0;
1088	break;
1089	}
1090
1091	// Iterate
1092	problemStatus_=-1;
1093	whileIterating(originalModel,sequenceIn,&info,crucialSj,phase);
1094	}
1095	}
1096	// clean up
1097	delete [] saveSolution;
1098	delete [] pair;
1099	finish();
1100	restoreData(data);
1101	return problemStatus_;
1102	}
1103	/*
1104	Reasons to come out:
1105	-1 iterations etc
1106	-2 inaccuracy
1107	-3 slight inaccuracy (and done iterations)
1108	-4 end of values pass and done iterations
1109	+0 looks optimal (might be infeasible - but we will investigate)
1110	+2 looks unbounded
1111	+3 max iterations
1112	*/
1113	int
1114	ClpSimplexPrimalQuadratic::whileIterating(
1115	const ClpSimplexPrimalQuadratic * originalModel,
1116	int & sequenceIn,
1117	ClpQuadraticInfo * info,
1118	int & crucialSj,
1119	int phase)
1120	{
1121	checkComplementary (info);
1122
1123	int returnCode=-1;
1124	double saveObjective = objectiveValue_;
1125	int numberXColumns = originalModel->numberColumns();
1126	int oldSequenceIn=sequenceIn;
1127	// status stays at -1 while iterating, >=0 finished, -2 to invert
1128	// status -3 to go to top without an invert
1129	while (problemStatus_==-1) {
1130	#ifdef CLP_DEBUG
1131	{
1132	int i;
1133	// not [1] as has information
1134	for (i=0;i<4;i++) {
1135	if (i!=1)
1136	rowArray_[i]->checkClear();
1137	}
1138	for (i=0;i<2;i++) {
1139	columnArray_[i]->checkClear();
1140	}
1141	}
1142	#endif
1143	// choose column to come in
1144	// can use pivotRow_ to update weights
1145	// pass in list of cost changes so can do row updates (rowArray_[1])
1146	// NOTE rowArray_[0] is used by computeDuals which is a
1147	// slow way of getting duals but might be used
1148	// Initially Dantzig and look at s variables
1149	// Only do if one not already chosen
1150	bool cleanupIteration;
1151	if (phase==2) {
1152	// values pass
1153	if (sequenceIn<0) {
1154	// get next
1155	int iColumn;
1156	int iStart = oldSequenceIn+1;
1157	for (iColumn=iStart;iColumn<numberXColumns;iColumn++) {
1158	double value = solution_[iColumn];
1159	double lower = lower_[iColumn];
1160	double upper = upper_[iColumn];
1161	if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
1162	if (getColumnStatus(iColumn)!=basic) {
1163	sequenceIn=iColumn;
1164	break;
1165	}
1166	}
1167	}
1168	if (sequenceIn<0) {
1169	iStart=max(iStart,numberColumns_);
1170	int numberXRows = originalModel->numberRows();
1171	for (iColumn=numberColumns_;iColumn<numberColumns_+numberXRows;
1172	iColumn++) {
1173	double value = solution_[iColumn];
1174	double lower = lower_[iColumn];
1175	double upper = upper_[iColumn];
1176	if (value>lower+primalTolerance_&&value<upper-primalTolerance_) {
1177	if (getColumnStatus(iColumn)!=basic) {
1178	sequenceIn=iColumn;
1179	break;
1180	}
1181	}
1182	}
1183	}
1184	}
1185	oldSequenceIn=sequenceIn;
1186	sequenceIn_ = sequenceIn;
1187	cleanupIteration=false;
1188	dualIn_ = solution_[sequenceIn_+numberRows_];
1189	valueIn_=solution_[sequenceIn_];
1190	if (dualIn_>0.0)
1191	directionIn_ = -1;
1192	else
1193	directionIn_ = 1;
1194	} else {
1195	if (sequenceIn<0) {
1196	primalColumn(rowArray_[1],rowArray_[2],rowArray_[3],
1197	columnArray_[0],columnArray_[1]);
1198	cleanupIteration=false;
1199	} else {
1200	sequenceIn_ = sequenceIn;
1201	cleanupIteration=true;
1202	}
1203	}
1204	pivotRow_=-1;
1205	sequenceOut_=-1;
1206	rowArray_[1]->clear();
1207	if (sequenceIn_>=0) {
1208	sequenceIn=-1;
1209	// we found a pivot column
1210	int chosen = sequenceIn_;
1211	// do second half of iteration
1212	while (chosen>=0) {
1213	objectiveValue_=saveObjective;
1214	returnCode=-1;
1215	pivotRow_=-1;
1216	sequenceOut_=-1;
1217	rowArray_[1]->clear();
1218	// we found a pivot column
1219	// update the incoming column
1220	sequenceIn_=chosen;
1221	chosen=-1;
1222	unpack(rowArray_[1]);
1223	factorization_->updateColumnFT(rowArray_[2],rowArray_[1]);
1224	if (cleanupIteration) {
1225	// move back to a version of primalColumn?
1226	valueIn_=solution_[sequenceIn_];
1227	// should keep pivot row of crucialSj as well (speed)
1228	int iSjRow=-1;
1229	{
1230	double * work=rowArray_[1]->denseVector();
1231	int number=rowArray_[1]->getNumElements();
1232	int * which=rowArray_[1]->getIndices();
1233	double tj = 0.0;
1234	int iIndex;
1235	for (iIndex=0;iIndex<number;iIndex++) {
1236	int iRow = which[iIndex];
1237	double alpha = work[iRow];
1238	int iPivot=pivotVariable_[iRow];
1239	if (iPivot==crucialSj) {
1240	tj = alpha;
1241	iSjRow = iRow;
1242	double d2 = solution_[crucialSj]/tj;
1243	// see which way to go
1244	if (d2>0)
1245	dj_[sequenceIn_]= -1.0;
1246	else
1247	dj_[sequenceIn_]= 1.0;
1248	break;
1249	}
1250	}
1251	if (!tj) {
1252	printf("trouble\n");
1253	assert (sequenceIn_>numberXColumns&&sequenceIn_<numberColumns_);
1254	dj_[sequenceIn_]=solution_[sequenceIn_];
1255	//assert(tj);
1256	}
1257	}
1258
1259	dualIn_=dj_[sequenceIn_];
1260	if (sequenceIn_>numberXColumns&&
1261	sequenceIn_<-numberColumns_) {
1262	// We can let flip to 0.0
1263	assert (cleanupIteration);
1264	if (solution_[sequenceIn_]>0.0) {
1265	nonLinearCost_->setBounds(sequenceIn_, 0.0,
1266	0.5*COIN_DBL_MAX);
1267	lower_[sequenceIn_]=0.0;
1268	upper_[sequenceIn_]=0.5*COIN_DBL_MAX;
1269	} else {
1270	nonLinearCost_->setBounds(sequenceIn_, -0.5*COIN_DBL_MAX,
1271	0.0);
1272	lower_[sequenceIn_]=-0.5*COIN_DBL_MAX;
1273	upper_[sequenceIn_]=0.0;
1274	}
1275	} else {
1276	// Let go through bounds
1277	nonLinearCost_->setBounds(sequenceIn_, -0.5*COIN_DBL_MAX,
1278	0.5*COIN_DBL_MAX);
1279	lower_[sequenceIn_]=-0.5*COIN_DBL_MAX;
1280	upper_[sequenceIn_]=0.5*COIN_DBL_MAX;
1281	}
1282	lowerIn_=lower_[sequenceIn_];
1283	upperIn_=upper_[sequenceIn_];
1284	if (dualIn_>0.0)
1285	directionIn_ = -1;
1286	else
1287	directionIn_ = 1;
1288	} else {
1289	if (sequenceIn_<numberColumns_) {
1290	// Set dj as value of slack
1291	crucialSj = sequenceIn_+numberRows_; // sj which should go to 0.0
1292	//dualIn_=solution_[crucialSj];
1293	} else {
1294	// Set dj as value of pi
1295	crucialSj = sequenceIn_-numberColumns_+numberXColumns; // pi which should go to 0.0
1296	//dualIn_=solution_[crucialSj];
1297	}
1298	}
1299	// save reduced cost
1300	//double saveDj = dualIn_;
1301	// do ratio test and re-compute dj
1302	// Note second parameter long enough for columns
1303	int result=primalRow(rowArray_[1],rowArray_[3],rowArray_[2],rowArray_[0],
1304	info,
1305	originalModel,crucialSj,cleanupIteration);
1306	saveObjective = objectiveValue_;
1307	if (pivotRow_>=0) {
1308	// if stable replace in basis
1309	int updateStatus = factorization_->replaceColumn(rowArray_[2],
1310	pivotRow_,
1311	alpha_);
1312	// if no pivots, bad update but reasonable alpha - take and invert
1313	if (updateStatus==2&&
1314	lastGoodIteration_==numberIterations_&&fabs(alpha_)>1.0e-5)
1315	updateStatus=4;
1316	if (updateStatus==1\|\|updateStatus==4) {
1317	// slight error
1318	if (factorization_->pivots()>5\|\|updateStatus==4) {
1319	returnCode=-3;
1320	}
1321	} else if (updateStatus==2) {
1322	// major error
1323	// better to have small tolerance even if slower
1324	factorization_->zeroTolerance(1.0e-15);
1325	int maxFactor = factorization_->maximumPivots();
1326	if (maxFactor>10) {
1327	if (forceFactorization_<0)
1328	forceFactorization_= maxFactor;
1329	forceFactorization_ = max (1,(forceFactorization_>>1));
1330	}
1331	// later we may need to unwind more e.g. fake bounds
1332	if(lastGoodIteration_ != numberIterations_) {
1333	rowArray_[1]->clear();
1334	pivotRow_=-1;
1335	returnCode=-4;
1336	// retry on return
1337	sequenceIn = sequenceIn_;
1338	break;
1339	} else {
1340	// need to reject something
1341	char x = isColumn(sequenceIn_) ? 'C' :'R';
1342	handler_->message(CLP_SIMPLEX_FLAG,messages_)
1343	<<x<<sequenceWithin(sequenceIn_)
1344	<<CoinMessageEol;
1345	setFlagged(sequenceIn_);
1346	lastBadIteration_ = numberIterations_; // say be more cautious
1347	rowArray_[1]->clear();
1348	pivotRow_=-1;
1349	returnCode = -5;
1350	break;
1351
1352	}
1353	} else if (updateStatus==3) {
1354	// out of memory
1355	// increase space if not many iterations
1356	if (factorization_->pivots()<
1357	0.5*factorization_->maximumPivots()&&
1358	factorization_->pivots()<200)
1359	factorization_->areaFactor(
1360	factorization_->areaFactor() * 1.1);
1361	returnCode =-2; // factorize now
1362	}
1363	// here do part of steepest - ready for next iteration
1364	primalColumnPivot_->updateWeights(rowArray_[1]);
1365	} else {
1366	if (pivotRow_==-1) {
1367	// no outgoing row is valid
1368	rowArray_[0]->clear();
1369	if (!factorization_->pivots()) {
1370	returnCode = 2; //say looks unbounded
1371	// do ray
1372	primalRay(rowArray_[1]);
1373	} else {
1374	returnCode = 4; //say looks unbounded but has iterated
1375	}
1376	break;
1377	} else {
1378	// flipping from bound to bound
1379	}
1380	}
1381
1382
1383	// update primal solution
1384
1385	double objectiveChange=0.0;
1386	// Cost on pivot row may change - may need to change dualIn
1387	double oldCost=0.0;
1388	if (pivotRow_>=0)
1389	oldCost = cost(pivotVariable_[pivotRow_]);
1390	// rowArray_[1] is not empty - used to update djs
1391	updatePrimalsInPrimal(rowArray_[1],theta_, objectiveChange);
1392	if (pivotRow_>=0)
1393	dualIn_ += (oldCost-cost(pivotVariable_[pivotRow_]));
1394	double oldValue = valueIn_;
1395	if (directionIn_==-1) {
1396	// as if from upper bound
1397	if (sequenceIn_!=sequenceOut_) {
1398	// variable becoming basic
1399	valueIn_ -= fabs(theta_);
1400	} else {
1401	valueIn_=lowerIn_;
1402	}
1403	} else {
1404	// as if from lower bound
1405	if (sequenceIn_!=sequenceOut_) {
1406	// variable becoming basic
1407	valueIn_ += fabs(theta_);
1408	} else {
1409	valueIn_=upperIn_;
1410	}
1411	}
1412	objectiveChange += dualIn_*(valueIn_-oldValue);
1413	// outgoing
1414	if (sequenceIn_!=sequenceOut_) {
1415	if (directionOut_>0) {
1416	valueOut_ = lowerOut_;
1417	} else {
1418	valueOut_ = upperOut_;
1419	}
1420	double lowerValue = lower_[sequenceOut_];
1421	double upperValue = upper_[sequenceOut_];
1422	assert(valueOut_>=lowerValue-primalTolerance_&&
1423	valueOut_<=upperValue+primalTolerance_);
1424	// may not be exactly at bound and bounds may have changed
1425	if (valueOut_<=lowerValue+primalTolerance_) {
1426	directionOut_=1;
1427	} else if (valueOut_>=upperValue-primalTolerance_) {
1428	directionOut_=-1;
1429	} else {
1430	printf("*** variable wandered off bound %g %g %g!\n",
1431	lowerValue,valueOut_,upperValue);
1432	if (valueOut_-lowerValue<=upperValue-valueOut_) {
1433	valueOut_=lowerValue;
1434	directionOut_=1;
1435	} else {
1436	valueOut_=upperValue;
1437	directionOut_=-1;
1438	}
1439	}
1440	solution_[sequenceOut_]=valueOut_;
1441	nonLinearCost_->setOne(sequenceOut_,valueOut_);
1442	}
1443	// change cost and bounds on incoming if primal
1444	nonLinearCost_->setOne(sequenceIn_,valueIn_);
1445	int whatNext=housekeeping(objectiveChange);
1446	checkComplementary (info);
1447
1448	if (whatNext==1) {
1449	returnCode =-2; // refactorize
1450	} else if (whatNext==2) {
1451	// maximum iterations or equivalent
1452	returnCode=3;
1453	} else if(numberIterations_ == lastGoodIteration_
1454	+ 2 * factorization_->maximumPivots()) {
1455	// done a lot of flips - be safe
1456	returnCode =-2; // refactorize
1457	}
1458	// may need to go round again
1459	cleanupIteration = true;
1460	double newDj;
1461	// may not be correct on second time
1462	if (sequenceIn_<numberXColumns)
1463	newDj = solution_[sequenceIn_+numberRows_];
1464	else
1465	newDj = solution_[sequenceIn_-numberColumns_+numberXColumns];
1466	newDj=1.0; // force
1467	if (result&&fabs(newDj)>dualTolerance_) {
1468	assert(sequenceOut_<numberXColumns\|\|
1469	sequenceOut_>=numberColumns_);
1470	if (sequenceOut_<numberXColumns) {
1471	chosen =sequenceOut_ + numberRows_; // sj variable
1472	} else {
1473	// Does this mean we can change pi
1474	int iRow = sequenceOut_-numberColumns_;
1475	if (iRow<numberRows_-numberXColumns) {
1476	int iPi = iRow+numberXColumns;
1477	printf("pi for row %d is %g\n",
1478	iRow,solution_[iPi]);
1479	chosen=iPi;
1480	} else {
1481	printf("Row %d is in column part\n",iRow);
1482	abort();
1483	}
1484	}
1485	} else {
1486	break;
1487	}
1488	}
1489	if (returnCode<-1&&returnCode>-5) {
1490	problemStatus_=-2; //
1491	} else if (returnCode==-5) {
1492	// something flagged - continue;
1493	} else if (returnCode==2) {
1494	problemStatus_=-5; // looks unbounded
1495	} else if (returnCode==4) {
1496	problemStatus_=-2; // looks unbounded but has iterated
1497	} else if (returnCode!=-1) {
1498	assert(returnCode==3);
1499	problemStatus_=3;
1500	}
1501	} else {
1502	// no pivot column
1503	#ifdef CLP_DEBUG
1504	if (handler_->logLevel()&32)
1505	printf("** no column pivot\n");
1506	#endif
1507	if (nonLinearCost_->numberInfeasibilities())
1508	problemStatus_=-4; // might be infeasible
1509	returnCode=0;
1510	break;
1511	}
1512	}
1513	return returnCode;
1514	}
1515	/*
1516	Row array has pivot column
1517	This chooses pivot row.
1518	For speed, we may need to go to a bucket approach when many
1519	variables go through bounds
1520	On exit rhsArray will have changes in costs of basic variables
1521	*/
1522	int
1523	ClpSimplexPrimalQuadratic::primalRow(CoinIndexedVector * rowArray,
1524	CoinIndexedVector * rhsArray,
1525	CoinIndexedVector * spareArray,
1526	CoinIndexedVector * spareArray2,
1527	ClpQuadraticInfo * info,
1528	const ClpSimplexPrimalQuadratic * originalModel,
1529	int crucialSj,
1530	bool cleanupIteration)
1531	{
1532	int result=1;
1533	// sequence stays as row number until end
1534	pivotRow_=-1;
1535	int numberSwapped=0;
1536	int numberRemaining=0;
1537
1538	int numberThru =0; // number gone thru a barrier
1539	int lastThru =0; // number gone thru a barrier on last time
1540
1541	double totalThru=0.0; // for when variables flip
1542	double acceptablePivot=1.0e-7;
1543	if (factorization_->pivots())
1544	acceptablePivot=1.0e-5; // if we have iterated be more strict
1545	double bestEverPivot=acceptablePivot;
1546	int lastPivotRow = -1;
1547	double lastPivot=0.0;
1548	double lastTheta=1.0e50;
1549	int lastNumberSwapped=0;
1550
1551	// use spareArrays to put ones looked at in
1552	// First one is list of candidates
1553	// We could compress if we really know we won't need any more
1554	// Second array has current set of pivot candidates
1555	// with a backup list saved in double * part of indexed vector
1556
1557	// for zeroing out arrays after
1558	int maximumSwapped=0;
1559	// pivot elements
1560	double * spare;
1561	// indices
1562	int * index, * indexSwapped;
1563	int * saveSwapped;
1564	spareArray->clear();
1565	spareArray2->clear();
1566	spare = spareArray->denseVector();
1567	index = spareArray->getIndices();
1568	saveSwapped = (int *) spareArray2->denseVector();
1569	indexSwapped = spareArray2->getIndices();
1570
1571	// we also need somewhere for effective rhs
1572	double * rhs=rhsArray->denseVector();
1573
1574	/*
1575	First we get a list of possible pivots. We can also see if the
1576	problem looks unbounded.
1577
1578	At first we increase theta and see what happens. We start
1579	theta at a reasonable guess. If in right area then we do bit by bit.
1580	We save possible pivot candidates
1581
1582	*/
1583
1584	// do first pass to get possibles
1585	// We can also see if unbounded
1586
1587	double * work=rowArray->denseVector();
1588	int number=rowArray->getNumElements();
1589	int * which=rowArray->getIndices();
1590
1591	// we need to swap sign if coming in from ub
1592	double way = directionIn_;
1593	double maximumMovement;
1594	if (way>0.0)
1595	maximumMovement = min(1.0e30,upperIn_-valueIn_);
1596	else
1597	maximumMovement = min(1.0e30,valueIn_-lowerIn_);
1598
1599	int iIndex;
1600
1601	// Work out coefficients for quadratic term
1602	// This is expanded one
1603	const CoinPackedMatrix * quadratic = quadraticObjective();
1604	const int * columnQuadratic = quadratic->getIndices();
1605	const int * columnQuadraticStart = quadratic->getVectorStarts();
1606	const int * columnQuadraticLength = quadratic->getVectorLengths();
1607	const double * quadraticElement = quadratic->getElements();
1608	const double * originalCost = originalModel->objective();
1609	// Use rhsArray
1610	rhsArray->clear();
1611	int * index2 = rhsArray->getIndices();
1612	int numberXColumns = originalModel->numberColumns();
1613	int number2=0;
1614	//int numberOriginalRows = originalModel->numberRows();
1615	// sj
1616	int iSjRow=-1;
1617	double tj = 0.0;
1618	for (iIndex=0;iIndex<number;iIndex++) {
1619	int iRow = which[iIndex];
1620	double alpha = -work[iRow]*way;
1621	int iPivot=pivotVariable_[iRow];
1622	if (iPivot<numberXColumns) {
1623	index2[number2++]=iPivot;
1624	rhs[iPivot]=alpha;
1625	//printf("col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
1626	} else {
1627	//printf("new col %d alpha %g solution %g\n",iPivot,alpha,solution_[iPivot]);
1628	if (iPivot==crucialSj) {
1629	tj = alpha;
1630	iSjRow = iRow;
1631	}
1632	}
1633	}
1634	// Incoming
1635	if (sequenceIn_<numberXColumns) {
1636	index2[number2++]=sequenceIn_;
1637	rhs[sequenceIn_]=way;
1638	printf("incoming col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
1639	} else {
1640	printf("incoming new col %d alpha %g solution %g\n",sequenceIn_,way,valueIn_);
1641	}
1642	rhsArray->setNumElements(number2);
1643	// Change in objective will be thetacoeff1 + thetatheta*coeff2
1644	double coeff1 = 0.0;
1645	double coeff2 = 0.0;
1646	if (numberIterations_>=0\|\|cleanupIteration) {
1647	for (iIndex=0;iIndex<number2;iIndex++) {
1648	int iColumn=index2[iIndex];
1649	//double valueI = solution_[iColumn];
1650	double alphaI = rhs[iColumn];
1651	coeff1 += alphaI*originalCost[iColumn];
1652	int j;
1653	for (j=columnQuadraticStart[iColumn];
1654	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1655	int jColumn = columnQuadratic[j];
1656	double valueJ = solution_[jColumn];
1657	double alphaJ = rhs[jColumn];
1658	double elementValue = quadraticElement[j];
1659	coeff1 += (valueJalphaI)elementValue;
1660	coeff2 += (alphaIalphaJ)elementValue;
1661	}
1662	}
1663	} else {
1664	const int * row = matrix_->getIndices();
1665	const int * columnStart = matrix_->getVectorStarts();
1666	const int * columnLength = matrix_->getVectorLengths();
1667	const double * element = matrix_->getElements();
1668	int j;
1669	int jRow = sequenceIn_+originalModel->numberRows();
1670	printf("sequence in %d, cost %g\n",sequenceIn_,
1671	upper_[jRow+numberColumns_]);
1672	double xx=0.0;
1673	for (j=0;j<numberColumns_;j++) {
1674	int k;
1675	for (k=columnStart[j];k<columnStart[j]+columnLength[j];k++) {
1676	if (row[k]==jRow) {
1677	printf ("col %d, el %g, sol %g, contr %g\n",
1678	j,element[k],solution_[j],element[k]*solution_[j]);
1679	xx+= element[k]*solution_[j];
1680	}
1681	}
1682	}
1683	printf("sum %g\n",xx);
1684	for (iIndex=0;iIndex<number2;iIndex++) {
1685	int iColumn=index2[iIndex];
1686	double valueI = solution_[iColumn];
1687	double alphaI = rhs[iColumn];
1688	//rhs[iColumn]=0.0;
1689	printf("Column %d, alpha %g sol %g cost %g, contr %g\n",
1690	iColumn,alphaI,valueI,originalCost[iColumn],
1691	alphaI*originalCost[iColumn]);
1692	coeff1 += alphaI*originalCost[iColumn];
1693	int j;
1694	for (j=columnQuadraticStart[iColumn];
1695	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
1696	int jColumn = columnQuadratic[j];
1697	double valueJ = solution_[jColumn];
1698	double alphaJ = rhs[jColumn];
1699	double elementValue = quadraticElement[j];
1700	printf("j %d alphaJ %g solJ %g el %g, contr %g\n",
1701	jColumn,alphaJ,valueJ,elementValue,
1702	(valueJalphaI)elementValue);
1703	coeff1 += (valueJalphaI)elementValue;
1704	coeff2 += (alphaIalphaJ)elementValue;
1705	}
1706	}
1707	}
1708	coeff2 *= 0.5;
1709	printf("coefficients %g %g - dualIn %g\n",coeff1,coeff2,dualIn_);
1710	if (!cleanupIteration)
1711	assert (fabs(way*coeff1-dualIn_)<1.0e-4);
1712	// interesting places are where derivative zero or sj goes to zero
1713	double d1,d2=1.0e50;
1714	if (fabs(coeff2)>1.0e-9)
1715	d1 = - 0.5*coeff1/coeff2;
1716	else if (coeff1<=1.0e-9)
1717	d1 = maximumMovement;
1718	else
1719	d1 = 0.0;
1720	if (fabs(tj)<1.0e-7) {
1721	if (sequenceIn_<numberXColumns)
1722	printf("column %d is basically linear\n",sequenceIn_);
1723	//assert(!columnQuadraticLength[sequenceIn_]);
1724	} else {
1725	d2 = -solution_[crucialSj]/tj;
1726	if (d2<0.0) {
1727	printf("d2 would be negative at %g\n",d2);
1728	d2=1.0e50;
1729	}
1730	}
1731	printf("derivative zero at %g, sj zero at %g\n",d1,d2);
1732	if (d1>1.0e10&&d2>1.0e10) {
1733	// linear variable entering
1734	// maybe we should have done dual iteration to force sj to 0.0
1735	printf("linear variable\n");
1736	}
1737	maximumMovement = min(maximumMovement,d1);
1738	maximumMovement = min(maximumMovement,d2);
1739	d2 = min(d1,d2);
1740
1741	rhsArray->clear();
1742	double tentativeTheta = maximumMovement;
1743	double upperTheta = maximumMovement;
1744
1745
1746	for (iIndex=0;iIndex<number;iIndex++) {
1747
1748	int iRow = which[iIndex];
1749	double alpha = work[iRow];
1750	int iPivot=pivotVariable_[iRow];
1751	alpha *= way;
1752	double oldValue = solution(iPivot);
1753	// get where in bound sequence
1754	if (alpha>0.0) {
1755	// basic variable going towards lower bound
1756	double bound = lower(iPivot);
1757	oldValue -= bound;
1758	} else if (alpha<0.0) {
1759	// basic variable going towards upper bound
1760	double bound = upper(iPivot);
1761	oldValue = bound-oldValue;
1762	}
1763	double value = oldValue-tentativeTheta*fabs(alpha);
1764	assert (oldValue>=-primalTolerance_*1.0001);
1765	if (value<-primalTolerance_) {
1766	// add to list
1767	spare[numberRemaining]=alpha;
1768	rhs[iRow]=oldValue;
1769	index[numberRemaining++]=iRow;
1770	double value=oldValue-upperTheta*fabs(alpha);
1771	if (value<-primalTolerance_&&fabs(alpha)>=acceptablePivot)
1772	upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1773	}
1774	}
1775
1776	// we need to keep where rhs non-zeros are
1777	int numberInRhs=numberRemaining;
1778	memcpy(rhsArray->getIndices(),index,numberInRhs*sizeof(int));
1779	rhsArray->setNumElements(numberInRhs);
1780
1781	theta_=maximumMovement;
1782
1783	bool goBackOne = false;
1784
1785	if (numberRemaining) {
1786
1787	// looks like pivoting
1788	// now try until reasonable theta
1789	tentativeTheta = max(10.0*upperTheta,1.0e-7);
1790	tentativeTheta = min(tentativeTheta,maximumMovement);
1791
1792	// loops increasing tentative theta until can't go through
1793
1794	while (tentativeTheta <= maximumMovement) {
1795	double thruThis = 0.0;
1796
1797	double bestPivot=acceptablePivot;
1798	pivotRow_ = -1;
1799
1800	numberSwapped = 0;
1801
1802	upperTheta = maximumMovement;
1803
1804	for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1805
1806	int iRow = index[iIndex];
1807	double alpha = spare[iIndex];
1808	double oldValue = rhs[iRow];
1809	double value = oldValue-tentativeTheta*fabs(alpha);
1810
1811	if (value<-primalTolerance_) {
1812	// how much would it cost to go thru
1813	thruThis += alpha*
1814	nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1815	// goes on swapped list (also means candidates if too many)
1816	indexSwapped[numberSwapped++]=iRow;
1817	if (fabs(alpha)>bestPivot) {
1818	bestPivot=fabs(alpha);
1819	pivotRow_ = iRow;
1820	theta_ = oldValue/bestPivot;
1821	}
1822	} else {
1823	value = oldValue-upperTheta*fabs(alpha);
1824	if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
1825	upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1826	}
1827	}
1828
1829	maximumSwapped = max(maximumSwapped,numberSwapped);
1830
1831	double dualCheck = - 2.0coeff2tentativeTheta - coeff1 - 99999999;
1832	// but make a bit more pessimistic
1833	dualCheck=max(dualCheck-100.0dualTolerance_,0.99dualCheck);
1834	if (totalThru+thruThis>=dualCheck) {
1835	// We should be pivoting in this batch
1836	// so compress down to this lot
1837
1838	int saveNumber = numberRemaining;
1839	numberRemaining=0;
1840	for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1841	int iRow = indexSwapped[iIndex];
1842	spare[numberRemaining]=way*work[iRow];
1843	index[numberRemaining++]=iRow;
1844	}
1845	memset(spare+numberRemaining,0,
1846	(saveNumber-numberRemaining)*sizeof(double));
1847	int iTry;
1848	#define MAXTRY 100
1849	// first get ratio with tolerance
1850	for (iTry=0;iTry<MAXTRY;iTry++) {
1851
1852	upperTheta=maximumMovement;
1853	numberSwapped = 0;
1854
1855	for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1856
1857	int iRow = index[iIndex];
1858	double alpha = fabs(spare[iIndex]);
1859	double oldValue = rhs[iRow];
1860	double value = oldValue-upperTheta*alpha;
1861
1862	if (value<-primalTolerance_ && alpha>=acceptablePivot)
1863	upperTheta = (oldValue+primalTolerance_)/alpha;
1864
1865	}
1866
1867	// now look at best in this lot
1868	bestPivot=acceptablePivot;
1869	pivotRow_=-1;
1870	for (iIndex=0;iIndex<numberRemaining;iIndex++) {
1871
1872	int iRow = index[iIndex];
1873	double alpha = spare[iIndex];
1874	double oldValue = rhs[iRow];
1875	double value = oldValue-upperTheta*fabs(alpha);
1876
1877	if (value<=0) {
1878	// how much would it cost to go thru
1879	totalThru += alpha*
1880	nonLinearCost_->changeInCost(pivotVariable_[iRow],alpha);
1881	// goes on swapped list (also means candidates if too many)
1882	indexSwapped[numberSwapped++]=iRow;
1883	if (fabs(alpha)>bestPivot) {
1884	bestPivot=fabs(alpha);
1885	theta_ = oldValue/bestPivot;
1886	pivotRow_=iRow;
1887	}
1888	} else {
1889	value = oldValue-upperTheta*fabs(alpha);
1890	if (value<-primalTolerance_ && fabs(alpha)>=acceptablePivot)
1891	upperTheta = (oldValue+primalTolerance_)/fabs(alpha);
1892	}
1893	}
1894
1895	maximumSwapped = max(maximumSwapped,numberSwapped);
1896	if (bestPivot<0.1*bestEverPivot&&
1897	bestEverPivot>1.0e-6&&bestPivot<1.0e-3) {
1898	// back to previous one
1899	goBackOne = true;
1900	break;
1901	} else if (pivotRow_==-1&&upperTheta>largeValue_) {
1902	if (lastPivot>acceptablePivot) {
1903	// back to previous one
1904	goBackOne = true;
1905	//break;
1906	} else {
1907	// can only get here if all pivots so far too small
1908	}
1909	break;
1910	} else {
1911	dualCheck = - 2.0coeff2theta_ - coeff1-9999999;
1912	if (totalThru>=dualCheck) {
1913	break; // no point trying another loop
1914	} else {
1915	// skip this lot
1916	nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1917	lastPivotRow=pivotRow_;
1918	lastTheta = theta_;
1919	lastThru = numberThru;
1920	numberThru += numberSwapped;
1921	lastNumberSwapped = numberSwapped;
1922	memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1923	if (bestPivot>bestEverPivot)
1924	bestEverPivot=bestPivot;
1925	}
1926	}
1927	}
1928	break;
1929	} else {
1930	// skip this lot
1931	nonLinearCost_->goThru(numberSwapped,way,indexSwapped, work,rhs);
1932	lastPivotRow=pivotRow_;
1933	lastTheta = theta_;
1934	lastThru = numberThru;
1935	numberThru += numberSwapped;
1936	lastNumberSwapped = numberSwapped;
1937	memcpy(saveSwapped,indexSwapped,lastNumberSwapped*sizeof(int));
1938	if (bestPivot>bestEverPivot)
1939	bestEverPivot=bestPivot;
1940	totalThru += thruThis;
1941	tentativeTheta = 2.0*upperTheta;
1942	}
1943	}
1944	// can get here without pivotRow_ set but with lastPivotRow
1945	if (goBackOne\|\|(pivotRow_<0&&lastPivotRow>=0)) {
1946	// back to previous one
1947	pivotRow_=lastPivotRow;
1948	theta_ = lastTheta;
1949	// undo this lot
1950	nonLinearCost_->goBack(lastNumberSwapped,saveSwapped,rhs);
1951	memcpy(indexSwapped,saveSwapped,lastNumberSwapped*sizeof(int));
1952	numberSwapped = lastNumberSwapped;
1953	}
1954	}
1955
1956	if (pivotRow_>=0) {
1957
1958	#define MINIMUMTHETA 1.0e-12
1959	// will we need to increase tolerance
1960	#ifdef CLP_DEBUG
1961	bool found=false;
1962	#endif
1963	double largestInfeasibility = primalTolerance_;
1964	if (theta_<MINIMUMTHETA) {
1965	theta_=MINIMUMTHETA;
1966	for (iIndex=0;iIndex<numberSwapped;iIndex++) {
1967	int iRow = indexSwapped[iIndex];
1968	#ifdef CLP_DEBUG
1969	if (iRow==pivotRow_)
1970	found=true;
1971	#endif
1972	largestInfeasibility = max (largestInfeasibility,
1973	-(rhs[iRow]-fabs(work[iRow])*theta_));
1974	}
1975	#ifdef CLP_DEBUG
1976	assert(found);
1977	if (largestInfeasibility>primalTolerance_&&(handler_->logLevel()&32))
1978	printf("Primal tolerance increased from %g to %g\n",
1979	primalTolerance_,largestInfeasibility);
1980	#endif
1981	primalTolerance_ = max(primalTolerance_,largestInfeasibility);
1982	}
1983	alpha_ = work[pivotRow_];
1984	// translate to sequence
1985	sequenceOut_ = pivotVariable_[pivotRow_];
1986	valueOut_ = solution(sequenceOut_);
1987	lowerOut_=lower_[sequenceOut_];
1988	upperOut_=upper_[sequenceOut_];
1989
1990	if (way<0.0)
1991	theta_ = - theta_;
1992	double newValue = valueOut_ - theta_*alpha_;
1993	if (alpha_*way<0.0) {
1994	directionOut_=-1; // to upper bound
1995	if (fabs(theta_)>1.0e-6)
1996	upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
1997	else
1998	upperOut_ = newValue;
1999	} else {
2000	directionOut_=1; // to lower bound
2001	if (fabs(theta_)>1.0e-6)
2002	lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2003	else
2004	lowerOut_ = newValue;
2005	}
2006	dualOut_ = reducedCost(sequenceOut_);
2007	} else {
2008	double trueMaximumMovement;
2009	if (way>0.0)
2010	trueMaximumMovement = min(1.0e30,upperIn_-valueIn_);
2011	else
2012	trueMaximumMovement = min(1.0e30,valueIn_-lowerIn_);
2013	if (maximumMovement<1.0e20&&maximumMovement==trueMaximumMovement) {
2014	// flip
2015	theta_ = maximumMovement;
2016	pivotRow_ = -2; // so we can tell its a flip
2017	result=0;
2018	sequenceOut_ = sequenceIn_;
2019	valueOut_ = valueIn_;
2020	dualOut_ = dualIn_;
2021	lowerOut_ = lowerIn_;
2022	upperOut_ = upperIn_;
2023	alpha_ = 0.0;
2024	if (way<0.0) {
2025	directionOut_=1; // to lower bound
2026	theta_ = lowerOut_ - valueOut_;
2027	} else {
2028	directionOut_=-1; // to upper bound
2029	theta_ = upperOut_ - valueOut_;
2030	}
2031	// we may still have sj to get rid of
2032	} else if (fabs(maximumMovement-d2)<dualTolerance_) {
2033	// sj going to zero
2034	result=0;
2035	assert (pivotRow_<0);
2036	nonLinearCost_->setBounds(crucialSj, 0.0,0.0);
2037	lower_[crucialSj]=0.0;
2038	upper_[crucialSj]=0.0;
2039	setColumnStatus(crucialSj,isFixed);
2040	pivotRow_ = iSjRow;
2041	alpha_ = work[pivotRow_];
2042	// translate to sequence
2043	sequenceOut_ = pivotVariable_[pivotRow_];
2044	valueOut_ = solution(sequenceOut_);
2045	lowerOut_=lower_[sequenceOut_];
2046	upperOut_=upper_[sequenceOut_];
2047	theta_ = d2;
2048	if (way<0.0)
2049	theta_ = - theta_;
2050	double newValue = valueOut_ - theta_*alpha_;
2051	if (alpha_*way<0.0) {
2052	directionOut_=-1; // to upper bound
2053	if (fabs(theta_)>1.0e-6)
2054	upperOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2055	else
2056	upperOut_ = newValue;
2057	} else {
2058	directionOut_=1; // to lower bound
2059	if (fabs(theta_)>1.0e-6)
2060	lowerOut_ = nonLinearCost_->nearest(sequenceOut_,newValue);
2061	else
2062	lowerOut_ = newValue;
2063	}
2064	//????
2065	dualOut_ = reducedCost(sequenceOut_);
2066	} else {
2067	// need to do something
2068	abort();
2069	}
2070	}
2071
2072	// clear arrays
2073
2074	memset(spare,0,numberRemaining*sizeof(double));
2075	memset(saveSwapped,0,maximumSwapped*sizeof(int));
2076
2077	// put back original bounds etc
2078	nonLinearCost_->goBackAll(rhsArray);
2079
2080	rhsArray->clear();
2081	// Change in objective will be thetacoeff1 + thetatheta*coeff2
2082	objectiveValue_ += theta_coeff1+theta_theta_*coeff2;
2083	printf("New objective value %g\n",objectiveValue_);
2084	{
2085	int iColumn;
2086	objectiveValue_ =0.0;
2087	CoinPackedMatrix * quadratic = originalModel->quadraticObjective();
2088	const int * columnQuadratic = quadratic->getIndices();
2089	const int * columnQuadraticStart = quadratic->getVectorStarts();
2090	const int * columnQuadraticLength = quadratic->getVectorLengths();
2091	const double * quadraticElement = quadratic->getElements();
2092	int numberColumns = originalModel->numberColumns();
2093	const double * objective = originalModel->objective();
2094	for (iColumn=0;iColumn<numberColumns;iColumn++)
2095	objectiveValue_ += objective[iColumn]*solution_[iColumn];
2096	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2097	double valueI = solution_[iColumn];
2098	int j;
2099	for (j=columnQuadraticStart[iColumn];
2100	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
2101	int jColumn = columnQuadratic[j];
2102	double valueJ = solution_[jColumn];
2103	double elementValue = quadraticElement[j];
2104	objectiveValue_ += 0.5valueIvalueJ*elementValue;
2105	}
2106	}
2107	printf("Objective value %g\n",objectiveValue());
2108	}
2109	return result;
2110	}
2111	// Just for debug
2112	int
2113	ClpSimplexPrimalQuadratic::checkComplementary (const
2114	ClpQuadraticInfo * info)
2115	{
2116	const ClpSimplex * originalModel= info->originalModel();
2117	int numberXColumns = originalModel->numberColumns();
2118	int i;
2119	for (i=0;i<numberXColumns;i++) {
2120	int jSequence = i+ numberRows_;
2121	if (getColumnStatus(i)==basic) {
2122	if (getColumnStatus(jSequence)==basic)
2123	printf("Struct %d (%g) and %d (%g) both basic\n",
2124	i,solution_[i],jSequence,solution_[jSequence]);
2125	}
2126	}
2127	int originalNumberRows = originalModel->numberRows();
2128	int offset = numberXColumns;
2129	for (i=0;i<originalNumberRows;i++) {
2130	int jSequence = i+offset;
2131	if (getRowStatus(i)==basic) {
2132	if (getColumnStatus(jSequence)==basic)
2133	printf("Row %d (%g) and %d (%g) both basic\n",
2134	i,solution_[i],jSequence,solution_[jSequence]);
2135	}
2136	}
2137	return 0;
2138	}
2139
2140	/* This creates the large version of QP and
2141	fills in quadratic information
2142	*/
2143	ClpSimplexPrimalQuadratic *
2144	ClpSimplexPrimalQuadratic::makeQuadratic(ClpQuadraticInfo & info)
2145	{
2146
2147	// Get list of non linear columns
2148	CoinPackedMatrix * quadratic = quadraticObjective();
2149	if (!quadratic\|\|!quadratic->getNumElements()) {
2150	// no quadratic part
2151	return NULL;
2152	}
2153
2154	int numberColumns = this->numberColumns();
2155	double * columnLower = this->columnLower();
2156	double * columnUpper = this->columnUpper();
2157	double * objective = this->objective();
2158	double * solution = this->primalColumnSolution();
2159	double * dj = this->dualColumnSolution();
2160	double * pi = this->dualRowSolution();
2161
2162	int numberRows = this->numberRows();
2163	double * rowLower = this->rowLower();
2164	double * rowUpper = this->rowUpper();
2165
2166	// and elements
2167	CoinPackedMatrix * matrix = this->matrix();
2168	const int * row = matrix->getIndices();
2169	const int * columnStart = matrix->getVectorStarts();
2170	const double * element = matrix->getElements();
2171	const int * columnLength = matrix->getVectorLengths();
2172
2173	int iColumn;
2174	const int * columnQuadratic = quadratic->getIndices();
2175	const int * columnQuadraticStart = quadratic->getVectorStarts();
2176	const int * columnQuadraticLength = quadratic->getVectorLengths();
2177	const double * quadraticElement = quadratic->getElements();
2178	#if 0
2179	// deliberate bad solution
2180	// Change to use phase
2181	//double * saveO = new double[numberColumns];
2182	//memcpy(saveO,objective,numberColumns*sizeof(double));
2183	//memset(objective,0,numberColumns*sizeof(double));
2184	tempMessage messageHandler(this);;
2185	passInMessageHandler(&messageHandler);
2186	factorization()->maximumPivots(1);
2187	primal();
2188	CoinMessageHandler handler2;
2189	passInMessageHandler(&handler2);
2190	factorization()->maximumPivots(100);
2191	setMaximumIterations(1000);
2192	#endif
2193	//memcpy(objective,saveO,numberColumns*sizeof(double));
2194	// Get a feasible solution
2195	//printf("For testing - deliberate bad solution\n");
2196	//columnUpper[0]=0.0;
2197	//columnLower[0]=0.0;
2198	//quadraticSLP(50,1.0e-4);
2199	//primal(1);
2200	//columnUpper[0]=COIN_DBL_MAX;
2201
2202	// Create larger problem
2203	// First workout how many rows extra
2204	info=ClpQuadraticInfo(this);
2205	int numberQuadratic = info.numberQuadraticColumns();
2206	int newNumberRows = numberRows+numberQuadratic;
2207	int newNumberColumns = numberColumns + numberRows + numberQuadratic;
2208	int numberElements = 2*matrix->getNumElements()
2209	+2*quadratic->getNumElements()
2210	+ numberQuadratic;
2211	double * elements2 = new double[numberElements];
2212	int * start2 = new int[newNumberColumns+1];
2213	int * row2 = new int[numberElements];
2214	double * objective2 = new double[newNumberColumns];
2215	double * columnLower2 = new double[newNumberColumns];
2216	double * columnUpper2 = new double[newNumberColumns];
2217	double * rowLower2 = new double[newNumberRows];
2218	double * rowUpper2 = new double[newNumberRows];
2219	const int * which = info.quadraticSequence();
2220	const int * back = info.backSequence();
2221	memcpy(rowLower2,rowLower,numberRows*sizeof(double));
2222	memcpy(rowUpper2,rowUpper,numberRows*sizeof(double));
2223	int iRow;
2224	for (iRow=0;iRow<numberQuadratic;iRow++) {
2225	double cost = objective[back[iRow]];
2226	rowLower2[iRow+numberRows]=cost;
2227	rowUpper2[iRow+numberRows]=cost;
2228	}
2229	memset(objective2,0,newNumberColumns*sizeof(double));
2230	// Get a row copy of quadratic objective in standard format
2231	CoinPackedMatrix copyQ;
2232	copyQ.reverseOrderedCopyOf(*quadratic);
2233	const int * columnQ = copyQ.getIndices();
2234	const CoinBigIndex * rowStartQ = copyQ.getVectorStarts();
2235	const int * rowLengthQ = copyQ.getVectorLengths();
2236	const double * elementByRowQ = copyQ.getElements();
2237	// Move solution across
2238	double * solution2 = new double[newNumberColumns];
2239	memset(solution2,0,newNumberColumns*sizeof(double));
2240	newNumberColumns=0;
2241	numberElements=0;
2242	start2[0]=0;
2243	// x
2244	memcpy(dj,objective,numberColumns*sizeof(double));
2245	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2246	// Original matrix
2247	columnLower2[iColumn]=columnLower[iColumn];
2248	columnUpper2[iColumn]=columnUpper[iColumn];
2249	solution2[iColumn]=solution[iColumn];
2250	int j;
2251	for (j=columnStart[iColumn];
2252	j<columnStart[iColumn]+columnLength[iColumn];
2253	j++) {
2254	elements2[numberElements]=element[j];
2255	row2[numberElements++]=row[j];
2256	}
2257	if (which[iColumn]>=0) {
2258	// Quadratic and modify djs
2259	for (j=columnQuadraticStart[iColumn];
2260	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
2261	j++) {
2262	int jColumn = columnQuadratic[j];
2263	double value = quadraticElement[j];
2264	if (iColumn!=jColumn)
2265	value *= 0.5;
2266	dj[jColumn] += solution[iColumn]*value;
2267	elements2[numberElements]=-value;
2268	row2[numberElements++]=which[jColumn]+numberRows;
2269	}
2270	for (j=rowStartQ[iColumn];
2271	j<rowStartQ[iColumn]+rowLengthQ[iColumn];
2272	j++) {
2273	int jColumn = columnQ[j];
2274	double value = elementByRowQ[j];
2275	if (iColumn!=jColumn) {
2276	value *= 0.5;
2277	dj[jColumn] += solution[iColumn]*value;
2278	elements2[numberElements]=-value;
2279	row2[numberElements++]=which[jColumn]+numberRows;
2280	}
2281	}
2282	}
2283	start2[iColumn+1]=numberElements;
2284	}
2285	newNumberColumns=numberColumns;
2286	// pi
2287	// Get a row copy in standard format
2288	CoinPackedMatrix copy;
2289	copy.reverseOrderedCopyOf(*(this->matrix()));
2290	// get matrix data pointers
2291	const int * column = copy.getIndices();
2292	const CoinBigIndex * rowStart = copy.getVectorStarts();
2293	const int * rowLength = copy.getVectorLengths();
2294	const double * elementByRow = copy.getElements();
2295	for (iRow=0;iRow<numberRows;iRow++) {
2296	solution2[newNumberColumns]=pi[iRow];
2297	double value = pi[iRow];
2298	columnLower2[newNumberColumns]=-COIN_DBL_MAX;
2299	columnUpper2[newNumberColumns]=COIN_DBL_MAX;
2300	int j;
2301	for (j=rowStart[iRow];
2302	j<rowStart[iRow]+rowLength[iRow];
2303	j++) {
2304	double elementValue=elementByRow[j];
2305	int jColumn = column[j];
2306	elements2[numberElements]=elementValue;
2307	row2[numberElements++]=jColumn+numberRows;
2308	dj[jColumn]-= value*elementValue;
2309	}
2310	newNumberColumns++;
2311	start2[newNumberColumns]=numberElements;
2312	}
2313	// djs
2314	for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
2315	columnLower2[newNumberColumns]=-COIN_DBL_MAX;
2316	columnUpper2[newNumberColumns]=COIN_DBL_MAX;
2317	solution2[newNumberColumns]=dj[iColumn];
2318	elements2[numberElements]=1.0;
2319	row2[numberElements++]=back[iColumn]+numberRows;
2320	newNumberColumns++;
2321	start2[newNumberColumns]=numberElements;
2322	}
2323	// Create model
2324	ClpSimplex * model2 = new ClpSimplex(*this);
2325	model2->resize(0,0);
2326	model2->loadProblem(newNumberColumns,newNumberRows,
2327	start2,row2, elements2,
2328	columnLower2,columnUpper2,
2329	objective2,
2330	rowLower2,rowUpper2);
2331	delete [] objective2;
2332	delete [] rowLower2;
2333	delete [] rowUpper2;
2334	delete [] columnLower2;
2335	delete [] columnUpper2;
2336	// Now create expanded quadratic objective for use in primalRow
2337	// Later pack down in some way
2338	start2[0]=0;
2339	numberElements=0;
2340	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2341	// Quadratic
2342	int j;
2343	for (j=columnQuadraticStart[iColumn];
2344	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];
2345	j++) {
2346	int jColumn = columnQuadratic[j];
2347	double value = quadraticElement[j];
2348	if (iColumn!=jColumn)
2349	value *= 0.5;
2350	elements2[numberElements]=value;
2351	row2[numberElements++]=jColumn;
2352	}
2353	for (j=rowStartQ[iColumn];
2354	j<rowStartQ[iColumn]+rowLengthQ[iColumn];
2355	j++) {
2356	int jColumn = columnQ[j];
2357	double value = elementByRowQ[j];
2358	if (iColumn!=jColumn) {
2359	value *= 0.5;
2360	elements2[numberElements]=value;
2361	row2[numberElements++]=jColumn;
2362	}
2363	}
2364	start2[iColumn+1]=numberElements;
2365	}
2366	// and pad
2367	for (;iColumn<newNumberColumns;iColumn++)
2368	start2[iColumn+1]=numberElements;
2369	// Load up objective
2370	model2->loadQuadraticObjective(newNumberColumns,start2,row2,elements2);
2371	delete [] start2;
2372	delete [] row2;
2373	delete [] elements2;
2374	model2->allSlackBasis();
2375	model2->scaling(false);
2376	model2->setLogLevel(this->logLevel());
2377	// Move solution across
2378	memcpy(model2->primalColumnSolution(),solution2,
2379	newNumberColumns*sizeof(double));
2380	columnLower2 = model2->columnLower();
2381	columnUpper2 = model2->columnUpper();
2382	delete [] solution2;
2383	solution2 = model2->primalColumnSolution();
2384	// Compute row activities and check feasible
2385	double * rowSolution2 = model2->primalRowSolution();
2386	memset(rowSolution2,0,newNumberRows*sizeof(double));
2387	model2->times(1.0,solution2,rowSolution2);
2388	rowLower2 = model2->rowLower();
2389	rowUpper2 = model2->rowUpper();
2390	#if 0
2391	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2392	Status xStatus = getColumnStatus(iColumn);
2393	bool isSuperBasic;
2394	int iS = iColumn+newNumberRows;
2395	double value = solution2[iS];
2396	if (fabs(value)>dualTolerance_)
2397	isSuperBasic=true;
2398	else
2399	isSuperBasic=false;
2400	// For moment take all x out of basis
2401	// Does not seem right
2402	isSuperBasic=true;
2403	model2->setColumnStatus(iColumn,xStatus);
2404	if (xStatus==basic) {
2405	if (!isSuperBasic) {
2406	model2->setRowStatus(numberRows+iColumn,basic);
2407	model2->setColumnStatus(iS,superBasic);
2408	} else {
2409	model2->setRowStatus(numberRows+iColumn,isFixed);
2410	model2->setColumnStatus(iS,basic);
2411	model2->setColumnStatus(iColumn,superBasic);
2412	}
2413	} else {
2414	model2->setRowStatus(numberRows+iColumn,isFixed);
2415	model2->setColumnStatus(iS,basic);
2416	}
2417	}
2418	for (iRow=0;iRow<numberRows;iRow++) {
2419	Status rowStatus = getRowStatus(iRow);
2420	model2->setRowStatus(iRow,rowStatus);
2421	if (rowStatus!=basic) {
2422	model2->setColumnStatus(iRow+numberColumns,basic); // make dual basic
2423	}
2424	assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
2425	assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
2426	}
2427	// why ?? - take duals out and adjust
2428	for (iRow=0;iRow<numberRows;iRow++) {
2429	model2->setRowStatus(iRow,basic);
2430	model2->setColumnStatus(iRow+numberColumns,superBasic);
2431	solution2[iRow+numberColumns]=0.0;
2432	}
2433	#else
2434	for (iRow=0;iRow<numberRows;iRow++) {
2435	assert (rowSolution2[iRow]>=rowLower2[iRow]-primalTolerance_);
2436	assert (rowSolution2[iRow]<=rowUpper2[iRow]+primalTolerance_);
2437	model2->setRowStatus(iRow,basic);
2438	model2->setColumnStatus(iRow+numberColumns,superBasic);
2439	solution2[iRow+numberColumns]=0.0;
2440	}
2441	for (iColumn=numberRows+numberColumns;iColumn<newNumberColumns;iColumn++) {
2442	model2->setColumnStatus(iColumn,basic);
2443	model2->setRowStatus(iColumn-numberColumns,isFixed);
2444	}
2445	#endif
2446	memset(rowSolution2,0,newNumberRows*sizeof(double));
2447	model2->times(1.0,solution2,rowSolution2);
2448	for (iColumn=0;iColumn<numberQuadratic;iColumn++) {
2449	int iS = back[iColumn]+newNumberRows;
2450	int iRow = iColumn+numberRows;
2451	double value = rowSolution2[iRow];
2452	if (value>rowUpper2[iRow]) {
2453	rowSolution2[iRow] = rowUpper2[iRow];
2454	solution2[iS]-=value-rowUpper2[iRow];
2455	} else {
2456	rowSolution2[iRow] = rowLower2[iRow];
2457	solution2[iS]-=value-rowLower2[iRow];
2458	}
2459	}
2460
2461
2462	// See if any s sub j have wrong sign and/or use djs from infeasibility objective
2463	double objectiveOffset;
2464	getDblParam(ClpObjOffset,objectiveOffset);
2465	double objValue = -objectiveOffset;
2466	for (iColumn=0;iColumn<numberColumns;iColumn++)
2467	objValue += objective[iColumn]*solution2[iColumn];
2468	for (iColumn=0;iColumn<numberColumns;iColumn++) {
2469	double valueI = solution2[iColumn];
2470	int j;
2471	for (j=columnQuadraticStart[iColumn];
2472	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
2473	int jColumn = columnQuadratic[j];
2474	double valueJ = solution2[jColumn];
2475	double elementValue = quadraticElement[j];
2476	objValue += 0.5valueIvalueJ*elementValue;
2477	}
2478	}
2479	printf("Objective value %g\n",objValue);
2480	for (iColumn=0;iColumn<newNumberColumns;iColumn++)
2481	printf("%d %g\n",iColumn,solution2[iColumn]);
2482	return (ClpSimplexPrimalQuadratic *) model2;
2483	}
2484
2485	// This moves solution back and deletes information
2486	int
2487	ClpSimplexPrimalQuadratic::endQuadratic(ClpSimplexPrimalQuadratic * quadraticModel,
2488	ClpQuadraticInfo & info)
2489	{
2490	memcpy(dualRowSolution(),quadraticModel->primalColumnSolution()+numberColumns_,numberRows_*sizeof(double));
2491	const double * solution2 = quadraticModel->primalColumnSolution();
2492	memcpy(primalColumnSolution(),solution2,numberColumns_*sizeof(double));
2493	memset(quadraticModel->primalRowSolution(),0,
2494	quadraticModel->numberRows()*sizeof(double));
2495	quadraticModel->times(1.0,quadraticModel->primalColumnSolution(),
2496	quadraticModel->primalRowSolution());
2497	memcpy(dualColumnSolution(),
2498	quadraticModel->primalColumnSolution()+numberRows_+numberColumns_,
2499	numberColumns_*sizeof(double));
2500
2501	int iColumn;
2502	double objectiveOffset;
2503	getDblParam(ClpObjOffset,objectiveOffset);
2504	double objValue = -objectiveOffset;
2505	double * objective = this->objective();
2506	for (iColumn=0;iColumn<numberColumns_;iColumn++)
2507	objValue += objective[iColumn]*solution2[iColumn];
2508	CoinPackedMatrix * quadratic = quadraticObjective();
2509	if (quadratic) {
2510	const int * columnQuadratic = quadratic->getIndices();
2511	const int * columnQuadraticStart = quadratic->getVectorStarts();
2512	const int * columnQuadraticLength = quadratic->getVectorLengths();
2513	const double * quadraticElement = quadratic->getElements();
2514	for (iColumn=0;iColumn<numberColumns_;iColumn++) {
2515	double valueI = solution2[iColumn];
2516	if (fabs(valueI)>1.0e-5) {
2517	int djColumn = iColumn+numberRows_+numberColumns_;
2518	assert(solution2[djColumn]<1.0e-7);
2519	}
2520	int j;
2521	for (j=columnQuadraticStart[iColumn];
2522	j<columnQuadraticStart[iColumn]+columnQuadraticLength[iColumn];j++) {
2523	int jColumn = columnQuadratic[j];
2524	double valueJ = solution2[jColumn];
2525	double elementValue = quadraticElement[j];
2526	objValue += 0.5valueIvalueJ*elementValue;
2527	}
2528	}
2529	objectiveValue_ = objValue + objectiveOffset;
2530	}
2531	printf("Objective value %g\n",objValue);
2532	return 0;
2533	}
2534
2535	/// Default constructor.
2536	ClpQuadraticInfo::ClpQuadraticInfo()
2537	: originalModel_(NULL),
2538	quadraticSequence_(NULL),
2539	backSequence_(NULL),
2540	crucialSj_(-1),
2541	numberXRows_(-1),
2542	numberXColumns_(-1),
2543	numberQuadraticColumns_(0)
2544	{
2545	}
2546	// Constructor from original model
2547	ClpQuadraticInfo::ClpQuadraticInfo(const ClpSimplex * model)
2548	: originalModel_(model),
2549	quadraticSequence_(NULL),
2550	backSequence_(NULL),
2551	crucialSj_(-1),
2552	numberXRows_(-1),
2553	numberXColumns_(-1),
2554	numberQuadraticColumns_(0)
2555	{
2556	if (originalModel_) {
2557	numberXRows_ = originalModel_->numberRows();
2558	numberXColumns_ = originalModel_->numberColumns();
2559	quadraticSequence_ = new int[numberXColumns_];
2560	backSequence_ = new int[numberXColumns_];
2561	int i;
2562	numberQuadraticColumns_=numberXColumns_;
2563	for (i=0;i<numberXColumns_;i++) {
2564	quadraticSequence_[i]=i;
2565	backSequence_[i]=i;
2566	}
2567	}
2568	}
2569	// Destructor
2570	ClpQuadraticInfo:: ~ClpQuadraticInfo()
2571	{
2572	delete [] quadraticSequence_;
2573	delete [] backSequence_;
2574	}
2575	// Copy
2576	ClpQuadraticInfo::ClpQuadraticInfo(const ClpQuadraticInfo& rhs)
2577	: originalModel_(rhs.originalModel_),
2578	quadraticSequence_(NULL),
2579	backSequence_(NULL),
2580	crucialSj_(rhs.crucialSj_),
2581	numberXRows_(rhs.numberXRows_),
2582	numberXColumns_(rhs.numberXColumns_),
2583	numberQuadraticColumns_(rhs.numberQuadraticColumns_)
2584	{
2585	if (numberXColumns_) {
2586	quadraticSequence_ = new int[numberXColumns_];
2587	memcpy(quadraticSequence_,rhs.quadraticSequence_,
2588	numberXColumns_*sizeof(int));
2589	backSequence_ = new int[numberXColumns_];
2590	memcpy(backSequence_,rhs.backSequence_,
2591	numberXColumns_*sizeof(int));
2592	}
2593	}
2594	// Assignment
2595	ClpQuadraticInfo &
2596	ClpQuadraticInfo::operator=(const ClpQuadraticInfo&rhs)
2597	{
2598	if (this != &rhs) {
2599	originalModel_ = rhs.originalModel_;
2600	delete [] quadraticSequence_;
2601	quadraticSequence_ = NULL;
2602	delete [] backSequence_;
2603	backSequence_ = NULL;
2604	crucialSj_ = rhs.crucialSj_;
2605	numberXRows_ = rhs.numberXRows_;
2606	numberXColumns_ = rhs.numberXColumns_;
2607	numberQuadraticColumns_=rhs.numberQuadraticColumns_;
2608	if (numberXColumns_) {
2609	quadraticSequence_ = new int[numberXColumns_];
2610	memcpy(quadraticSequence_,rhs.quadraticSequence_,
2611	numberXColumns_*sizeof(int));
2612	backSequence_ = new int[numberXColumns_];
2613	memcpy(backSequence_,rhs.backSequence_,
2614	numberXColumns_*sizeof(int));
2615	}
2616	}
2617	return *this;
2618	}

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

Download in other formats: