ClpSimplexNonlinear.cpp source code [OGDF/src/coin/Clp/ClpSimplexNonlinear.cpp]

1	/ $Id: ClpSimplexNonlinear.cpp 1665 2011-01-04 17:55:54Z lou $ /
2	// Copyright (C) 2004, International Business Machines
3	// Corporation and others. All Rights Reserved.
4	// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6	#include "CoinPragma.hpp"
7
8	#include <math.h>
9	#include "CoinHelperFunctions.hpp"
10	#include "ClpHelperFunctions.hpp"
11	#include "ClpSimplexNonlinear.hpp"
12	#include "ClpFactorization.hpp"
13	#include "ClpNonLinearCost.hpp"
14	#include "ClpLinearObjective.hpp"
15	#include "ClpConstraint.hpp"
16	#include "ClpQuadraticObjective.hpp"
17	#include "CoinPackedMatrix.hpp"
18	#include "CoinIndexedVector.hpp"
19	#include "ClpPrimalColumnPivot.hpp"
20	#include "ClpMessage.hpp"
21	#include "ClpEventHandler.hpp"
22	#include <cfloat>
23	#include <cassert>
24	#include <string>
25	#include <stdio.h>
26	#include <iostream>
27	#ifndef NDEBUG
28	#define CLP_DEBUG 1
29	#endif
30	// primal
31	int ClpSimplexNonlinear::primal ()
32	{
33
34	int ifValuesPass = `1`;
35	algorithm_ = +`3`;
36
37	// save data
38	ClpDataSave data = saveData();
39	matrix_->refresh(this); // make sure matrix okay
40
41	// Save objective
42	ClpObjective * saveObjective = NULL;
43	if (objective_->type() > `1`) {
44	// expand to full if quadratic
45	#ifndef NO_RTTI
46	ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
47	#else
48	ClpQuadraticObjective * quadraticObj = NULL;
49	if (objective_->type() == `2`)
50	quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
51	#endif
52	// for moment only if no scaling
53	// May be faster if switched off - but can't see why
54	if (!quadraticObj->fullMatrix() && (!rowScale_ && !scalingFlag_) && objectiveScale_ == `1.0`) {
55	saveObjective = objective_;
56	objective_ = new ClpQuadraticObjective (*quadraticObj, `1`);
57	}
58	}
59	double bestObjectiveWhenFlagged = COIN_DBL_MAX;
60	int pivotMode = `15`;
61	//pivotMode=20;
62
63	// initialize - maybe values pass and algorithm_ is +1
64	// true does something (? perturbs)
65	if (!startup(true)) {
66
67	// Set average theta
68	nonLinearCost_->setAverageTheta(`1.0e3`);
69	int lastCleaned = `0`; // last time objective or bounds cleaned up
70
71	// Say no pivot has occurred (for steepest edge and updates)
72	pivotRow_ = -`2`;
73
74	// This says whether to restore things etc
75	int factorType = `0`;
76	// Start check for cycles
77	progress_.startCheck();
78	/*
79	Status of problem:
80	0 - optimal
81	1 - infeasible
82	2 - unbounded
83	-1 - iterating
84	-2 - factorization wanted
85	-3 - redo checking without factorization
86	-4 - looks infeasible
87	-5 - looks unbounded
88	*/
89	while (problemStatus_ < `0`) {
90	int iRow, iColumn;
91	// clear
92	for (iRow = `0`; iRow < `4`; iRow++) {
93	rowArray_[iRow]->clear();
94	}
95
96	for (iColumn = `0`; iColumn < `2`; iColumn++) {
97	columnArray_[iColumn]->clear();
98	}
99
100	// give matrix (and model costs and bounds a chance to be
101	// refreshed (normally null)
102	matrix_->refresh(this);
103	// If getting nowhere - why not give it a kick
104	// If we have done no iterations - special
105	if (lastGoodIteration_ == numberIterations_ && factorType)
106	factorType = `3`;
107
108	// may factorize, checks if problem finished
109	if (objective_->type() > `1` && lastFlaggedIteration_ >= `0` &&
110	numberIterations_ > lastFlaggedIteration_ + `507`) {
111	unflag();
112	lastFlaggedIteration_ = numberIterations_;
113	if (pivotMode >= `10`) {
114	pivotMode--;
115	#ifdef CLP_DEBUG
116	if (handler_->logLevel() & `32`)
117	printf("pivot mode now %d\n", pivotMode);
118	#endif
119	if (pivotMode == `9`)
120	pivotMode = `0`; // switch off fast attempt
121	}
122	}
123	statusOfProblemInPrimal(lastCleaned, factorType, &progress_, true,
124	bestObjectiveWhenFlagged);
125
126	// Say good factorization
127	factorType = `1`;
128
129	// Say no pivot has occurred (for steepest edge and updates)
130	pivotRow_ = -`2`;
131
132	// exit if victory declared
133	if (problemStatus_ >= `0`)
134	break;
135
136	// test for maximum iterations
137	if (hitMaximumIterations() \|\| (ifValuesPass == `2` && firstFree_ < `0`)) {
138	problemStatus_ = `3`;
139	break;
140	}
141
142	if (firstFree_ < `0`) {
143	if (ifValuesPass) {
144	// end of values pass
145	ifValuesPass = `0`;
146	int status = eventHandler_->event(ClpEventHandler::endOfValuesPass);
147	if (status >= `0`) {
148	problemStatus_ = `5`;
149	secondaryStatus_ = ClpEventHandler::endOfValuesPass;
150	break;
151	}
152	}
153	}
154	// Check event
155	{
156	int status = eventHandler_->event(ClpEventHandler::endOfFactorization);
157	if (status >= `0`) {
158	problemStatus_ = `5`;
159	secondaryStatus_ = ClpEventHandler::endOfFactorization;
160	break;
161	}
162	}
163	// Iterate
164	whileIterating(pivotMode);
165	}
166	}
167	// if infeasible get real values
168	if (problemStatus_ == `1`) {
169	infeasibilityCost_ = `0.0`;
170	createRim(`1` + `4`);
171	nonLinearCost_->checkInfeasibilities(`0.0`);
172	sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
173	numberPrimalInfeasibilities_ = nonLinearCost_->numberInfeasibilities();
174	// and get good feasible duals
175	computeDuals(NULL);
176	}
177	// correct objective value
178	if (numberColumns_)
179	objectiveValue_ = nonLinearCost_->feasibleCost() + objective_->nonlinearOffset();
180	objectiveValue_ /= (objectiveScale_ * rhsScale_);
181	// clean up
182	unflag();
183	finish();
184	restoreData(data);
185	// restore objective if full
186	if (saveObjective) {
187	delete objective_;
188	objective_ = saveObjective;
189	}
190	return problemStatus_;
191	}
192	/ Refactorizes if necessary*
193	Checks if finished. Updates status.
194	lastCleaned refers to iteration at which some objective/feasibility
195	cleaning too place.
196
197	type - 0 initial so set up save arrays etc
198	- 1 normal -if good update save
199	- 2 restoring from saved
200	*/
201	void
202	ClpSimplexNonlinear::statusOfProblemInPrimal(int & lastCleaned, int type,
203	ClpSimplexProgress * progress,
204	bool doFactorization,
205	double & bestObjectiveWhenFlagged)
206	{
207	int dummy; // for use in generalExpanded
208	if (type == `2`) {
209	// trouble - restore solution
210	CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_ );
211	CoinMemcpyN(savedSolution_ + numberColumns_ , numberRows_, rowActivityWork_);
212	CoinMemcpyN(savedSolution_ , numberColumns_, columnActivityWork_);
213	// restore extra stuff
214	matrix_->generalExpanded(this, `6`, dummy);
215	forceFactorization_ = `1`; // a bit drastic but ..
216	pivotRow_ = -`1`; // say no weights update
217	changeMade_++; // say change made
218	}
219	int saveThreshold = factorization_->sparseThreshold();
220	int tentativeStatus = problemStatus_;
221	int numberThrownOut = `1`; // to loop round on bad factorization in values pass
222	while (numberThrownOut) {
223	if (problemStatus_ > -`3` \|\| problemStatus_ == -`4`) {
224	// factorize
225	// later on we will need to recover from singularities
226	// also we could skip if first time
227	// do weights
228	// This may save pivotRow_ for use
229	if (doFactorization)
230	primalColumnPivot_->saveWeights(this, `1`);
231
232	if (type && doFactorization) {
233	// is factorization okay?
234	int factorStatus = internalFactorize(`1`);
235	if (factorStatus) {
236	if (type != `1` \|\| largestPrimalError_ > `1.0e3`
237	\|\| largestDualError_ > `1.0e3`) {
238	// was \|\|largestDualError_>1.0e3\|\|objective_->type()>1) {
239	// switch off dense
240	int saveDense = factorization_->denseThreshold();
241	factorization_->setDenseThreshold(`0`);
242	// make sure will do safe factorization
243	pivotVariable_[`0`] = -`1`;
244	internalFactorize(`2`);
245	factorization_->setDenseThreshold(saveDense);
246	// Go to safe
247	factorization_->pivotTolerance(`0.99`);
248	// restore extra stuff
249	matrix_->generalExpanded(this, `6`, dummy);
250	} else {
251	// no - restore previous basis
252	CoinMemcpyN(saveStatus_, (numberColumns_ + numberRows_), status_ );
253	CoinMemcpyN(savedSolution_ + numberColumns_ , numberRows_, rowActivityWork_);
254	CoinMemcpyN(savedSolution_ , numberColumns_, columnActivityWork_);
255	// restore extra stuff
256	matrix_->generalExpanded(this, `6`, dummy);
257	matrix_->generalExpanded(this, `5`, dummy);
258	forceFactorization_ = `1`; // a bit drastic but ..
259	type = `2`;
260	// Go to safe
261	factorization_->pivotTolerance(`0.99`);
262	if (internalFactorize(`1`) != `0`)
263	largestPrimalError_ = `1.0e4`; // force other type
264	}
265	// flag incoming
266	if (sequenceIn_ >= `0` && getStatus(sequenceIn_) != basic) {
267	setFlagged(sequenceIn_);
268	saveStatus_[sequenceIn_] = status_[sequenceIn_];
269	}
270	changeMade_++; // say change made
271	}
272	}
273	if (problemStatus_ != -`4`)
274	problemStatus_ = -`3`;
275	}
276	// at this stage status is -3 or -5 if looks unbounded
277	// get primal and dual solutions
278	// put back original costs and then check
279	createRim(`4`);
280	// May need to do more if column generation
281	dummy = `4`;
282	matrix_->generalExpanded(this, `9`, dummy);
283	numberThrownOut = gutsOfSolution(NULL, NULL, (firstFree_ >= `0`));
284	if (numberThrownOut) {
285	problemStatus_ = tentativeStatus;
286	doFactorization = true;
287	}
288	}
289	// Double check reduced costs if no action
290	if (progress->lastIterationNumber(`0`) == numberIterations_) {
291	if (primalColumnPivot_->looksOptimal()) {
292	numberDualInfeasibilities_ = `0`;
293	sumDualInfeasibilities_ = `0.0`;
294	}
295	}
296	// Check if looping
297	int loop;
298	if (type != `2`)
299	loop = progress->looping();
300	else
301	loop = -`1`;
302	if (loop >= `0`) {
303	if (!problemStatus_) {
304	// declaring victory
305	numberPrimalInfeasibilities_ = `0`;
306	sumPrimalInfeasibilities_ = `0.0`;
307	} else {
308	problemStatus_ = loop; //exit if in loop
309	problemStatus_ = `10`; // instead - try other algorithm
310	}
311	problemStatus_ = `10`; // instead - try other algorithm
312	return ;
313	} else if (loop < -`1`) {
314	// Is it time for drastic measures
315	if (nonLinearCost_->numberInfeasibilities() && progress->badTimes() > `5` &&
316	progress->oddState() < `10` && progress->oddState() >= `0`) {
317	progress->newOddState();
318	nonLinearCost_->zapCosts();
319	}
320	// something may have changed
321	gutsOfSolution(NULL, NULL, true);
322	}
323	// If progress then reset costs
324	if (loop == -`1` && !nonLinearCost_->numberInfeasibilities() && progress->oddState() < `0`) {
325	createRim(`4`, false); // costs back
326	delete nonLinearCost_;
327	nonLinearCost_ = new ClpNonLinearCost (this);
328	progress->endOddState();
329	gutsOfSolution(NULL, NULL, true);
330	}
331	// Flag to say whether to go to dual to clean up
332	bool goToDual = false;
333	// really for free variables in
334	//if((progressFlag_&2)!=0)
335	//problemStatus_=-1;
336	progressFlag_ = `0`; //reset progress flag
337
338	handler_->message(CLP_SIMPLEX_STATUS, messages_)
339	<< numberIterations_ << nonLinearCost_->feasibleReportCost();
340	handler_->printing(nonLinearCost_->numberInfeasibilities() > `0`)
341	<< nonLinearCost_->sumInfeasibilities() << nonLinearCost_->numberInfeasibilities();
342	handler_->printing(sumDualInfeasibilities_ > `0.0`)
343	<< sumDualInfeasibilities_ << numberDualInfeasibilities_;
344	handler_->printing(numberDualInfeasibilitiesWithoutFree_
345	< numberDualInfeasibilities_)
346	<< numberDualInfeasibilitiesWithoutFree_;
347	handler_->message() << CoinMessageEol;
348	if (!primalFeasible()) {
349	nonLinearCost_->checkInfeasibilities(primalTolerance_);
350	gutsOfSolution(NULL, NULL, true);
351	nonLinearCost_->checkInfeasibilities(primalTolerance_);
352	}
353	double trueInfeasibility = nonLinearCost_->sumInfeasibilities();
354	if (trueInfeasibility > `1.0`) {
355	// If infeasibility going up may change weights
356	double testValue = trueInfeasibility - `1.0e-4` * (`10.0` + trueInfeasibility);
357	if(progress->lastInfeasibility() < testValue) {
358	if (infeasibilityCost_ < `1.0e14`) {
359	infeasibilityCost_ *= `1.5`;
360	if (handler_->logLevel() == `63`)
361	printf("increasing weight to %g\n", infeasibilityCost_);
362	gutsOfSolution(NULL, NULL, true);
363	}
364	}
365	}
366	// we may wish to say it is optimal even if infeasible
367	bool alwaysOptimal = (specialOptions_ & `1`) != `0`;
368	// give code benefit of doubt
369	if (sumOfRelaxedDualInfeasibilities_ == `0.0` &&
370	sumOfRelaxedPrimalInfeasibilities_ == `0.0`) {
371	// say optimal (with these bounds etc)
372	numberDualInfeasibilities_ = `0`;
373	sumDualInfeasibilities_ = `0.0`;
374	numberPrimalInfeasibilities_ = `0`;
375	sumPrimalInfeasibilities_ = `0.0`;
376	}
377	// had \|\|(type==3&&problemStatus_!=-5) -- ??? why ????
378	if (dualFeasible() \|\| problemStatus_ == -`4`) {
379	// see if extra helps
380	if (nonLinearCost_->numberInfeasibilities() &&
381	(nonLinearCost_->sumInfeasibilities() > `1.0e-3` \|\| sumOfRelaxedPrimalInfeasibilities_)
382	&& !alwaysOptimal) {
383	//may need infeasiblity cost changed
384	// we can see if we can construct a ray
385	// make up a new objective
386	double saveWeight = infeasibilityCost_;
387	// save nonlinear cost as we are going to switch off costs
388	ClpNonLinearCost * nonLinear = nonLinearCost_;
389	// do twice to make sure Primal solution has settled
390	// put non-basics to bounds in case tolerance moved
391	// put back original costs
392	createRim(`4`);
393	nonLinearCost_->checkInfeasibilities(`0.0`);
394	gutsOfSolution(NULL, NULL, true);
395
396	infeasibilityCost_ = `1.0e100`;
397	// put back original costs
398	createRim(`4`);
399	nonLinearCost_->checkInfeasibilities(primalTolerance_);
400	// may have fixed infeasibilities - double check
401	if (nonLinearCost_->numberInfeasibilities() == `0`) {
402	// carry on
403	problemStatus_ = -`1`;
404	infeasibilityCost_ = saveWeight;
405	nonLinearCost_->checkInfeasibilities(primalTolerance_);
406	} else {
407	nonLinearCost_ = NULL;
408	// scale
409	int i;
410	for (i = `0`; i < numberRows_ + numberColumns_; i++)
411	cost_[i] *= `1.0e-95`;
412	gutsOfSolution(NULL, NULL, false);
413	nonLinearCost_ = nonLinear;
414	infeasibilityCost_ = saveWeight;
415	if ((infeasibilityCost_ >= `1.0e18` \|\|
416	numberDualInfeasibilities_ == `0`) && perturbation_ == `101`) {
417	goToDual = unPerturb(); // stop any further perturbation
418	if (nonLinearCost_->sumInfeasibilities() > `1.0e-1`)
419	goToDual = false;
420	nonLinearCost_->checkInfeasibilities(primalTolerance_);
421	numberDualInfeasibilities_ = `1`; // carry on
422	problemStatus_ = -`1`;
423	}
424	if (infeasibilityCost_ >= `1.0e20` \|\|
425	numberDualInfeasibilities_ == `0`) {
426	// we are infeasible - use as ray
427	delete [] ray_;
428	ray_ = new double [numberRows_];
429	CoinMemcpyN(dual_, numberRows_, ray_);
430	// and get feasible duals
431	infeasibilityCost_ = `0.0`;
432	createRim(`4`);
433	nonLinearCost_->checkInfeasibilities(primalTolerance_);
434	gutsOfSolution(NULL, NULL, true);
435	// so will exit
436	infeasibilityCost_ = `1.0e30`;
437	// reset infeasibilities
438	sumPrimalInfeasibilities_ = nonLinearCost_->sumInfeasibilities();
439	numberPrimalInfeasibilities_ =
440	nonLinearCost_->numberInfeasibilities();
441	}
442	if (infeasibilityCost_ < `1.0e20`) {
443	infeasibilityCost_ *= `5.0`;
444	changeMade_++; // say change made
445	unflag();
446	handler_->message(CLP_PRIMAL_WEIGHT, messages_)
447	<< infeasibilityCost_
448	<< CoinMessageEol;
449	// put back original costs and then check
450	createRim(`4`);
451	nonLinearCost_->checkInfeasibilities(`0.0`);
452	gutsOfSolution(NULL, NULL, true);
453	problemStatus_ = -`1`; //continue
454	goToDual = false;
455	} else {
456	// say infeasible
457	problemStatus_ = `1`;
458	}
459	}
460	} else {
461	// may be optimal
462	if (perturbation_ == `101`) {
463	goToDual = unPerturb(); // stop any further perturbation
464	lastCleaned = -`1`; // carry on
465	}
466	bool unflagged = unflag() != `0`;
467	if ( lastCleaned != numberIterations_ \|\| unflagged) {
468	handler_->message(CLP_PRIMAL_OPTIMAL, messages_)
469	<< primalTolerance_
470	<< CoinMessageEol;
471	double current = nonLinearCost_->feasibleReportCost();
472	if (numberTimesOptimal_ < `4`) {
473	if (bestObjectiveWhenFlagged <= current) {
474	numberTimesOptimal_++;
475	#ifdef CLP_DEBUG
476	if (handler_->logLevel() & `32`)
477	printf("%d times optimal, current %g best %g\n", numberTimesOptimal_,
478	current, bestObjectiveWhenFlagged);
479	#endif
480	} else {
481	bestObjectiveWhenFlagged = current;
482	}
483	changeMade_++; // say change made
484	if (numberTimesOptimal_ == `1`) {
485	// better to have small tolerance even if slower
486	factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), `1.0e-15`));
487	}
488	lastCleaned = numberIterations_;
489	if (primalTolerance_ != dblParam_[ClpPrimalTolerance])
490	handler_->message(CLP_PRIMAL_ORIGINAL, messages_)
491	<< CoinMessageEol;
492	double oldTolerance = primalTolerance_;
493	primalTolerance_ = dblParam_[ClpPrimalTolerance];
494	// put back original costs and then check
495	createRim(`4`);
496	nonLinearCost_->checkInfeasibilities(oldTolerance);
497	gutsOfSolution(NULL, NULL, true);
498	if (sumOfRelaxedDualInfeasibilities_ == `0.0` &&
499	sumOfRelaxedPrimalInfeasibilities_ == `0.0`) {
500	// say optimal (with these bounds etc)
501	numberDualInfeasibilities_ = `0`;
502	sumDualInfeasibilities_ = `0.0`;
503	numberPrimalInfeasibilities_ = `0`;
504	sumPrimalInfeasibilities_ = `0.0`;
505	}
506	if (dualFeasible() && !nonLinearCost_->numberInfeasibilities() && lastCleaned >= `0`)
507	problemStatus_ = `0`;
508	else
509	problemStatus_ = -`1`;
510	} else {
511	problemStatus_ = `0`; // optimal
512	if (lastCleaned < numberIterations_) {
513	handler_->message(CLP_SIMPLEX_GIVINGUP, messages_)
514	<< CoinMessageEol;
515	}
516	}
517	} else {
518	problemStatus_ = `0`; // optimal
519	}
520	}
521	} else {
522	// see if looks unbounded
523	if (problemStatus_ == -`5`) {
524	if (nonLinearCost_->numberInfeasibilities()) {
525	if (infeasibilityCost_ > `1.0e18` && perturbation_ == `101`) {
526	// back off weight
527	infeasibilityCost_ = `1.0e13`;
528	unPerturb(); // stop any further perturbation
529	}
530	//we need infeasiblity cost changed
531	if (infeasibilityCost_ < `1.0e20`) {
532	infeasibilityCost_ *= `5.0`;
533	changeMade_++; // say change made
534	unflag();
535	handler_->message(CLP_PRIMAL_WEIGHT, messages_)
536	<< infeasibilityCost_
537	<< CoinMessageEol;
538	// put back original costs and then check
539	createRim(`4`);
540	gutsOfSolution(NULL, NULL, true);
541	problemStatus_ = -`1`; //continue
542	} else {
543	// say unbounded
544	problemStatus_ = `2`;
545	}
546	} else {
547	// say unbounded
548	problemStatus_ = `2`;
549	}
550	} else {
551	if(type == `3` && problemStatus_ != -`5`)
552	unflag(); // odd
553	// carry on
554	problemStatus_ = -`1`;
555	}
556	}
557	// save extra stuff
558	matrix_->generalExpanded(this, `5`, dummy);
559	if (type == `0` \|\| type == `1`) {
560	if (type != `1` \|\| !saveStatus_) {
561	// create save arrays
562	delete [] saveStatus_;
563	delete [] savedSolution_;
564	saveStatus_ = new unsigned char [numberRows_+numberColumns_];
565	savedSolution_ = new double [numberRows_+numberColumns_];
566	}
567	// save arrays
568	CoinMemcpyN(status_, (numberColumns_ + numberRows_), saveStatus_);
569	CoinMemcpyN(rowActivityWork_, numberRows_, savedSolution_ + numberColumns_ );
570	CoinMemcpyN(columnActivityWork_, numberColumns_, savedSolution_ );
571	}
572	if (doFactorization) {
573	// restore weights (if saved) - also recompute infeasibility list
574	if (tentativeStatus > -`3`)
575	primalColumnPivot_->saveWeights(this, (type < `2`) ? `2` : `4`);
576	else
577	primalColumnPivot_->saveWeights(this, `3`);
578	if (saveThreshold) {
579	// use default at present
580	factorization_->sparseThreshold(`0`);
581	factorization_->goSparse();
582	}
583	}
584	if (problemStatus_ < `0` && !changeMade_) {
585	problemStatus_ = `4`; // unknown
586	}
587	lastGoodIteration_ = numberIterations_;
588	//if (goToDual)
589	//problemStatus_=10; // try dual
590	// Allow matrices to be sorted etc
591	int fake = -`999`; // signal sort
592	matrix_->correctSequence(this, fake, fake);
593	}
594	/*
595	Reasons to come out:
596	-1 iterations etc
597	-2 inaccuracy
598	-3 slight inaccuracy (and done iterations)
599	-4 end of values pass and done iterations
600	+0 looks optimal (might be infeasible - but we will investigate)
601	+2 looks unbounded
602	+3 max iterations
603	*/
604	int
605	ClpSimplexNonlinear::whileIterating(int & pivotMode)
606	{
607	// Say if values pass
608	//int ifValuesPass=(firstFree_>=0) ? 1 : 0;
609	int ifValuesPass = `1`;
610	int returnCode = -`1`;
611	// status stays at -1 while iterating, >=0 finished, -2 to invert
612	// status -3 to go to top without an invert
613	int numberInterior = `0`;
614	int nextUnflag = `10`;
615	int nextUnflagIteration = numberIterations_ + `10`;
616	// get two arrays
617	double * array1 = new double[`2`*(numberRows_+numberColumns_)];
618	double solutionError = -`1.0`;
619	while (problemStatus_ == -`1`) {
620	int result;
621	rowArray_[`1`]->clear();
622	//#define CLP_DEBUG
623	#if CLP_DEBUG > 1
624	rowArray_[`0`]->checkClear();
625	rowArray_[`1`]->checkClear();
626	rowArray_[`2`]->checkClear();
627	rowArray_[`3`]->checkClear();
628	columnArray_[`0`]->checkClear();
629	#endif
630	if (numberInterior >= `5`) {
631	// this can go when we have better minimization
632	if (pivotMode < `10`)
633	pivotMode = `1`;
634	unflag();
635	#ifdef CLP_DEBUG
636	if (handler_->logLevel() & `32`)
637	printf("interior unflag\n");
638	#endif
639	numberInterior = `0`;
640	nextUnflag = `10`;
641	nextUnflagIteration = numberIterations_ + `10`;
642	} else {
643	if (numberInterior > nextUnflag &&
644	numberIterations_ > nextUnflagIteration) {
645	nextUnflagIteration = numberIterations_ + `10`;
646	nextUnflag += `10`;
647	unflag();
648	#ifdef CLP_DEBUG
649	if (handler_->logLevel() & `32`)
650	printf("unflagging as interior\n");
651	#endif
652	}
653	}
654	pivotRow_ = -`1`;
655	result = pivotColumn(rowArray_[`3`], rowArray_[`0`],
656	columnArray_[`0`], rowArray_[`1`], pivotMode, solutionError,
657	array1);
658	if (result) {
659	if (result == `2` && sequenceIn_ < `0`) {
660	// does not look good
661	double currentObj;
662	double thetaObj;
663	double predictedObj;
664	objective_->stepLength(this, solution_, solution_, `0.0`,
665	currentObj, thetaObj, predictedObj);
666	if (currentObj == predictedObj) {
667	#ifdef CLP_INVESTIGATE
668	printf("looks bad - no change in obj %g\n", currentObj);
669	#endif
670	if (factorization_->pivots())
671	result = `3`;
672	else
673	problemStatus_ = `0`;
674	}
675	}
676	if (result == `3`)
677	break; // null vector not accurate
678	#ifdef CLP_DEBUG
679	if (handler_->logLevel() & `32`) {
680	double currentObj;
681	double thetaObj;
682	double predictedObj;
683	objective_->stepLength(this, solution_, solution_, `0.0`,
684	currentObj, thetaObj, predictedObj);
685	printf("obj %g after interior move\n", currentObj);
686	}
687	#endif
688	// just move and try again
689	if (pivotMode < `10`) {
690	pivotMode = result - `1`;
691	numberInterior++;
692	}
693	continue;
694	} else {
695	if (pivotMode < `10`) {
696	if (theta_ > `0.001`)
697	pivotMode = `0`;
698	else if (pivotMode == `2`)
699	pivotMode = `1`;
700	}
701	numberInterior = `0`;
702	nextUnflag = `10`;
703	nextUnflagIteration = numberIterations_ + `10`;
704	}
705	sequenceOut_ = -`1`;
706	rowArray_[`1`]->clear();
707	if (sequenceIn_ >= `0`) {
708	// we found a pivot column
709	assert (!flagged(sequenceIn_));
710	#ifdef CLP_DEBUG
711	if ((handler_->logLevel() & `32`)) {
712	char x = isColumn(sequenceIn_) ? `'C'` : `'R'`;
713	std::cout << "pivot column " <<
714	x << sequenceWithin(sequenceIn_) << std::endl;
715	}
716	#endif
717	// do second half of iteration
718	if (pivotRow_ < `0` && theta_ < `1.0e-8`) {
719	assert (sequenceIn_ >= `0`);
720	returnCode = pivotResult(ifValuesPass);
721	} else {
722	// specialized code
723	returnCode = pivotNonlinearResult();
724	//printf("odd pivrow %d\n",sequenceOut_);
725	if (sequenceOut_ >= `0` && theta_ < `1.0e-5`) {
726	if (getStatus(sequenceOut_) != isFixed) {
727	if (getStatus(sequenceOut_) == atUpperBound)
728	solution_[sequenceOut_] = upper_[sequenceOut_];
729	else if (getStatus(sequenceOut_) == atLowerBound)
730	solution_[sequenceOut_] = lower_[sequenceOut_];
731	setFlagged(sequenceOut_);
732	}
733	}
734	}
735	if (returnCode < -`1` && returnCode > -`5`) {
736	problemStatus_ = -`2`; //
737	} else if (returnCode == -`5`) {
738	// something flagged - continue;
739	} else if (returnCode == `2`) {
740	problemStatus_ = -`5`; // looks unbounded
741	} else if (returnCode == `4`) {
742	problemStatus_ = -`2`; // looks unbounded but has iterated
743	} else if (returnCode != -`1`) {
744	assert(returnCode == `3`);
745	problemStatus_ = `3`;
746	}
747	} else {
748	// no pivot column
749	#ifdef CLP_DEBUG
750	if (handler_->logLevel() & `32`)
751	printf("** no column pivot\n");
752	#endif
753	if (pivotMode < `10`) {
754	// looks optimal
755	primalColumnPivot_->setLooksOptimal(true);
756	} else {
757	pivotMode--;
758	#ifdef CLP_DEBUG
759	if (handler_->logLevel() & `32`)
760	printf("pivot mode now %d\n", pivotMode);
761	#endif
762	if (pivotMode == `9`)
763	pivotMode = `0`; // switch off fast attempt
764	unflag();
765	}
766	if (nonLinearCost_->numberInfeasibilities())
767	problemStatus_ = -`4`; // might be infeasible
768	returnCode = `0`;
769	break;
770	}
771	}
772	delete [] array1;
773	return returnCode;
774	}
775	// Creates direction vector
776	void
777	ClpSimplexNonlinear::directionVector (CoinIndexedVector * vectorArray,
778	CoinIndexedVector * spare1, CoinIndexedVector * spare2,
779	int pivotMode2,
780	double & normFlagged, double & normUnflagged,
781	int & numberNonBasic)
782	{
783	#if CLP_DEBUG > 1
784	vectorArray->checkClear();
785	spare1->checkClear();
786	spare2->checkClear();
787	#endif
788	double *array = vectorArray->denseVector();
789	int * index = vectorArray->getIndices();
790	int number = `0`;
791	sequenceIn_ = -`1`;
792	normFlagged = `0.0`;
793	normUnflagged = `1.0`;
794	double dualTolerance2 = CoinMin(`1.0e-8`, `1.0e-2` * dualTolerance_);
795	double dualTolerance3 = CoinMin(`1.0e-2`, `1.0e3` * dualTolerance_);
796	if (!numberNonBasic) {
797	//if (nonLinearCost_->sumInfeasibilities()>1.0e-4)
798	//printf("infeasible\n");
799	if (!pivotMode2 \|\| pivotMode2 >= `10`) {
800	normUnflagged = `0.0`;
801	double bestDj = `0.0`;
802	double bestSuper = `0.0`;
803	double sumSuper = `0.0`;
804	sequenceIn_ = -`1`;
805	int nSuper = `0`;
806	for (int iSequence = `0`; iSequence < numberColumns_ + numberRows_; iSequence++) {
807	array[iSequence] = `0.0`;
808	if (flagged(iSequence)) {
809	// accumulate norm
810	switch(getStatus(iSequence)) {
811
812	case basic:
813	case ClpSimplex::isFixed:
814	break;
815	case atUpperBound:
816	if (dj_[iSequence] > dualTolerance3)
817	normFlagged += dj_[iSequence] * dj_[iSequence];
818	break;
819	case atLowerBound:
820	if (dj_[iSequence] < -dualTolerance3)
821	normFlagged += dj_[iSequence] * dj_[iSequence];
822	break;
823	case isFree:
824	case superBasic:
825	if (fabs(dj_[iSequence]) > dualTolerance3)
826	normFlagged += dj_[iSequence] * dj_[iSequence];
827	break;
828	}
829	continue;
830	}
831	switch(getStatus(iSequence)) {
832
833	case basic:
834	case ClpSimplex::isFixed:
835	break;
836	case atUpperBound:
837	if (dj_[iSequence] > dualTolerance_) {
838	if (dj_[iSequence] > dualTolerance3)
839	normUnflagged += dj_[iSequence] * dj_[iSequence];
840	if (pivotMode2 < `10`) {
841	array[iSequence] = -dj_[iSequence];
842	index[number++] = iSequence;
843	} else {
844	if (dj_[iSequence] > bestDj) {
845	bestDj = dj_[iSequence];
846	sequenceIn_ = iSequence;
847	}
848	}
849	}
850	break;
851	case atLowerBound:
852	if (dj_[iSequence] < -dualTolerance_) {
853	if (dj_[iSequence] < -dualTolerance3)
854	normUnflagged += dj_[iSequence] * dj_[iSequence];
855	if (pivotMode2 < `10`) {
856	array[iSequence] = -dj_[iSequence];
857	index[number++] = iSequence;
858	} else {
859	if (-dj_[iSequence] > bestDj) {
860	bestDj = -dj_[iSequence];
861	sequenceIn_ = iSequence;
862	}
863	}
864	}
865	break;
866	case isFree:
867	case superBasic:
868	//#define ALLSUPER
869	#define NOSUPER
870	#ifndef ALLSUPER
871	if (fabs(dj_[iSequence]) > dualTolerance_) {
872	if (fabs(dj_[iSequence]) > dualTolerance3)
873	normUnflagged += dj_[iSequence] * dj_[iSequence];
874	nSuper++;
875	bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper);
876	sumSuper += fabs(dj_[iSequence]);
877	}
878	if (fabs(dj_[iSequence]) > dualTolerance2) {
879	array[iSequence] = -dj_[iSequence];
880	index[number++] = iSequence;
881	}
882	#else
883	array[iSequence] = -dj_[iSequence];
884	index[number++] = iSequence;
885	if (pivotMode2 >= `10`)
886	bestSuper = CoinMax(fabs(dj_[iSequence]), bestSuper);
887	#endif
888	break;
889	}
890	}
891	#ifdef NOSUPER
892	// redo
893	bestSuper = sumSuper;
894	if(sequenceIn_ >= `0` && bestDj > bestSuper) {
895	int j;
896	// get rid of superbasics
897	for (j = `0`; j < number; j++) {
898	int iSequence = index[j];
899	array[iSequence] = `0.0`;
900	}
901	number = `0`;
902	array[sequenceIn_] = -dj_[sequenceIn_];
903	index[number++] = sequenceIn_;
904	} else {
905	sequenceIn_ = -`1`;
906	}
907	#else
908	if (pivotMode2 >= `10` \|\| !nSuper) {
909	bool takeBest = true;
910	if (pivotMode2 == `100` && nSuper > `1`)
911	takeBest = false;
912	if(sequenceIn_ >= `0` && takeBest) {
913	if (fabs(dj_[sequenceIn_]) > bestSuper) {
914	array[sequenceIn_] = -dj_[sequenceIn_];
915	index[number++] = sequenceIn_;
916	} else {
917	sequenceIn_ = -`1`;
918	}
919	} else {
920	sequenceIn_ = -`1`;
921	}
922	}
923	#endif
924	#ifdef CLP_DEBUG
925	if (handler_->logLevel() & `32`) {
926	if (sequenceIn_ >= `0`)
927	printf("%d superBasic, chosen %d - dj %g\n", nSuper, sequenceIn_,
928	dj_[sequenceIn_]);
929	else
930	printf("%d superBasic - none chosen\n", nSuper);
931	}
932	#endif
933	} else {
934	double bestDj = `0.0`;
935	double saveDj = `0.0`;
936	if (sequenceOut_ >= `0`) {
937	saveDj = dj_[sequenceOut_];
938	dj_[sequenceOut_] = `0.0`;
939	switch(getStatus(sequenceOut_)) {
940
941	case basic:
942	sequenceOut_ = -`1`;
943	case ClpSimplex::isFixed:
944	break;
945	case atUpperBound:
946	if (dj_[sequenceOut_] > dualTolerance_) {
947	#ifdef CLP_DEBUG
948	if (handler_->logLevel() & `32`)
949	printf("after pivot out %d values %g %g %g, dj %g\n",
950	sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
951	upper_[sequenceOut_], dj_[sequenceOut_]);
952	#endif
953	}
954	break;
955	case atLowerBound:
956	if (dj_[sequenceOut_] < -dualTolerance_) {
957	#ifdef CLP_DEBUG
958	if (handler_->logLevel() & `32`)
959	printf("after pivot out %d values %g %g %g, dj %g\n",
960	sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
961	upper_[sequenceOut_], dj_[sequenceOut_]);
962	#endif
963	}
964	break;
965	case isFree:
966	case superBasic:
967	if (dj_[sequenceOut_] > dualTolerance_) {
968	#ifdef CLP_DEBUG
969	if (handler_->logLevel() & `32`)
970	printf("after pivot out %d values %g %g %g, dj %g\n",
971	sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
972	upper_[sequenceOut_], dj_[sequenceOut_]);
973	#endif
974	} else if (dj_[sequenceOut_] < -dualTolerance_) {
975	#ifdef CLP_DEBUG
976	if (handler_->logLevel() & `32`)
977	printf("after pivot out %d values %g %g %g, dj %g\n",
978	sequenceOut_, lower_[sequenceOut_], solution_[sequenceOut_],
979	upper_[sequenceOut_], dj_[sequenceOut_]);
980	#endif
981	}
982	break;
983	}
984	}
985	// Go for dj
986	pivotMode2 = `3`;
987	for (int iSequence = `0`; iSequence < numberColumns_ + numberRows_; iSequence++) {
988	array[iSequence] = `0.0`;
989	if (flagged(iSequence))
990	continue;
991	switch(getStatus(iSequence)) {
992
993	case basic:
994	case ClpSimplex::isFixed:
995	break;
996	case atUpperBound:
997	if (dj_[iSequence] > dualTolerance_) {
998	double distance = CoinMin(`1.0e-2`, solution_[iSequence] - lower_[iSequence]);
999	double merit = distance * dj_[iSequence];
1000	if (pivotMode2 == `1`)
1001	merit = `1.0e-20`; // discourage*
1002	if (pivotMode2 == `3`)
1003	merit = fabs(dj_[iSequence]);
1004	if (merit > bestDj) {
1005	sequenceIn_ = iSequence;
1006	bestDj = merit;
1007	}
1008	}
1009	break;
1010	case atLowerBound:
1011	if (dj_[iSequence] < -dualTolerance_) {
1012	double distance = CoinMin(`1.0e-2`, upper_[iSequence] - solution_[iSequence]);
1013	double merit = -distance * dj_[iSequence];
1014	if (pivotMode2 == `1`)
1015	merit = `1.0e-20`; // discourage*
1016	if (pivotMode2 == `3`)
1017	merit = fabs(dj_[iSequence]);
1018	if (merit > bestDj) {
1019	sequenceIn_ = iSequence;
1020	bestDj = merit;
1021	}
1022	}
1023	break;
1024	case isFree:
1025	case superBasic:
1026	if (dj_[iSequence] > dualTolerance_) {
1027	double distance = CoinMin(`1.0e-2`, solution_[iSequence] - lower_[iSequence]);
1028	distance = CoinMin(solution_[iSequence] - lower_[iSequence],
1029	upper_[iSequence] - solution_[iSequence]);
1030	double merit = distance * dj_[iSequence];
1031	if (pivotMode2 == `1`)
1032	merit = distance;
1033	if (pivotMode2 == `3`)
1034	merit = fabs(dj_[iSequence]);
1035	if (merit > bestDj) {
1036	sequenceIn_ = iSequence;
1037	bestDj = merit;
1038	}
1039	} else if (dj_[iSequence] < -dualTolerance_) {
1040	double distance = CoinMin(`1.0e-2`, upper_[iSequence] - solution_[iSequence]);
1041	distance = CoinMin(solution_[iSequence] - lower_[iSequence],
1042	upper_[iSequence] - solution_[iSequence]);
1043	double merit = -distance * dj_[iSequence];
1044	if (pivotMode2 == `1`)
1045	merit = distance;
1046	if (pivotMode2 == `3`)
1047	merit = fabs(dj_[iSequence]);
1048	if (merit > bestDj) {
1049	sequenceIn_ = iSequence;
1050	bestDj = merit;
1051	}
1052	}
1053	break;
1054	}
1055	}
1056	if (sequenceOut_ >= `0`) {
1057	dj_[sequenceOut_] = saveDj;
1058	sequenceOut_ = -`1`;
1059	}
1060	if (sequenceIn_ >= `0`) {
1061	array[sequenceIn_] = -dj_[sequenceIn_];
1062	index[number++] = sequenceIn_;
1063	}
1064	}
1065	numberNonBasic = number;
1066	} else {
1067	// compute norms
1068	normUnflagged = `0.0`;
1069	for (int iSequence = `0`; iSequence < numberColumns_ + numberRows_; iSequence++) {
1070	if (flagged(iSequence)) {
1071	// accumulate norm
1072	switch(getStatus(iSequence)) {
1073
1074	case basic:
1075	case ClpSimplex::isFixed:
1076	break;
1077	case atUpperBound:
1078	if (dj_[iSequence] > dualTolerance_)
1079	normFlagged += dj_[iSequence] * dj_[iSequence];
1080	break;
1081	case atLowerBound:
1082	if (dj_[iSequence] < -dualTolerance_)
1083	normFlagged += dj_[iSequence] * dj_[iSequence];
1084	break;
1085	case isFree:
1086	case superBasic:
1087	if (fabs(dj_[iSequence]) > dualTolerance_)
1088	normFlagged += dj_[iSequence] * dj_[iSequence];
1089	break;
1090	}
1091	}
1092	}
1093	// re-use list
1094	number = `0`;
1095	int j;
1096	for (j = `0`; j < numberNonBasic; j++) {
1097	int iSequence = index[j];
1098	if (flagged(iSequence))
1099	continue;
1100	switch(getStatus(iSequence)) {
1101
1102	case basic:
1103	case ClpSimplex::isFixed:
1104	continue; //abort();
1105	break;
1106	case atUpperBound:
1107	if (dj_[iSequence] > dualTolerance_) {
1108	number++;
1109	normUnflagged += dj_[iSequence] * dj_[iSequence];
1110	}
1111	break;
1112	case atLowerBound:
1113	if (dj_[iSequence] < -dualTolerance_) {
1114	number++;
1115	normUnflagged += dj_[iSequence] * dj_[iSequence];
1116	}
1117	break;
1118	case isFree:
1119	case superBasic:
1120	if (fabs(dj_[iSequence]) > dualTolerance_) {
1121	number++;
1122	normUnflagged += dj_[iSequence] * dj_[iSequence];
1123	}
1124	break;
1125	}
1126	array[iSequence] = -dj_[iSequence];
1127	}
1128	// switch to large
1129	normUnflagged = `1.0`;
1130	if (!number) {
1131	for (j = `0`; j < numberNonBasic; j++) {
1132	int iSequence = index[j];
1133	array[iSequence] = `0.0`;
1134	}
1135	numberNonBasic = `0`;
1136	}
1137	number = numberNonBasic;
1138	}
1139	if (number) {
1140	int j;
1141	// Now do basic ones - how do I compensate for small basic infeasibilities?
1142	int iRow;
1143	for (iRow = `0`; iRow < numberRows_; iRow++) {
1144	int iPivot = pivotVariable_[iRow];
1145	double value = `0.0`;
1146	if (solution_[iPivot] > upper_[iPivot]) {
1147	value = upper_[iPivot] - solution_[iPivot];
1148	} else if (solution_[iPivot] < lower_[iPivot]) {
1149	value = lower_[iPivot] - solution_[iPivot];
1150	}
1151	//if (value)
1152	//printf("inf %d %g %g %g\n",iPivot,lower_[iPivot],solution_[iPivot],
1153	// upper_[iPivot]);
1154	//value=0.0;
1155	value *= -`1.0e0`;
1156	if (value) {
1157	array[iPivot] = value;
1158	index[number++] = iPivot;
1159	}
1160	}
1161	double * array2 = spare1->denseVector();
1162	int * index2 = spare1->getIndices();
1163	int number2 = `0`;
1164	times(-`1.0`, array, array2);
1165	array = array + numberColumns_;
1166	for (iRow = `0`; iRow < numberRows_; iRow++) {
1167	double value = array2[iRow] + array[iRow];
1168	if (value) {
1169	array2[iRow] = value;
1170	index2[number2++] = iRow;
1171	} else {
1172	array2[iRow] = `0.0`;
1173	}
1174	}
1175	array -= numberColumns_;
1176	spare1->setNumElements(number2);
1177	// Ftran
1178	factorization_->updateColumn(spare2, spare1);
1179	number2 = spare1->getNumElements();
1180	for (j = `0`; j < number2; j++) {
1181	int iSequence = index2[j];
1182	double value = array2[iSequence];
1183	array2[iSequence] = `0.0`;
1184	if (value) {
1185	int iPivot = pivotVariable_[iSequence];
1186	double oldValue = array[iPivot];
1187	if (!oldValue) {
1188	array[iPivot] = value;
1189	index[number++] = iPivot;
1190	} else {
1191	// something already there
1192	array[iPivot] = value + oldValue;
1193	}
1194	}
1195	}
1196	spare1->setNumElements(`0`);
1197	}
1198	vectorArray->setNumElements(number);
1199	}
1200	#define MINTYPE 1
1201	#if MINTYPE==2
1202	static double
1203	innerProductIndexed(const double * region1, int size, const double * region2, const int * which)
1204	{
1205	int i;
1206	double value = `0.0`;
1207	for (i = `0`; i < size; i++) {
1208	int j = which[i];
1209	value += region1[j] * region2[j];
1210	}
1211	return value;
1212	}
1213	#endif
1214	/*
1215	Row array and column array have direction
1216	Returns 0 - can do normal iteration (basis change)
1217	1 - no basis change
1218	*/
1219	int
1220	ClpSimplexNonlinear::pivotColumn(CoinIndexedVector * longArray,
1221	CoinIndexedVector * rowArray,
1222	CoinIndexedVector * columnArray,
1223	CoinIndexedVector * spare,
1224	int & pivotMode,
1225	double & solutionError,
1226	double * dArray)
1227	{
1228	// say not optimal
1229	primalColumnPivot_->setLooksOptimal(false);
1230	double acceptablePivot = `1.0e-10`;
1231	int lastSequenceIn = -`1`;
1232	if (pivotMode && pivotMode < `10`) {
1233	acceptablePivot = `1.0e-6`;
1234	if (factorization_->pivots())
1235	acceptablePivot = `1.0e-5`; // if we have iterated be more strict
1236	}
1237	double acceptableBasic = `1.0e-7`;
1238
1239	int number = longArray->getNumElements();
1240	int numberTotal = numberRows_ + numberColumns_;
1241	int bestSequence = -`1`;
1242	int bestBasicSequence = -`1`;
1243	double eps = `1.0e-1`;
1244	eps = `1.0e-6`;
1245	double basicTheta = `1.0e30`;
1246	double objTheta = `0.0`;
1247	bool finished = false;
1248	sequenceIn_ = -`1`;
1249	int nPasses = `0`;
1250	int nTotalPasses = `0`;
1251	int nBigPasses = `0`;
1252	double djNorm0 = `0.0`;
1253	double djNorm = `0.0`;
1254	double normFlagged = `0.0`;
1255	double normUnflagged = `0.0`;
1256	int localPivotMode = pivotMode;
1257	bool allFinished = false;
1258	bool justOne = false;
1259	int returnCode = `1`;
1260	double currentObj;
1261	double predictedObj;
1262	double thetaObj;
1263	objective_->stepLength(this, solution_, solution_, `0.0`,
1264	currentObj, predictedObj, thetaObj);
1265	double saveObj = currentObj;
1266	#if MINTYPE ==2
1267	// try Shanno's method
1268	//would be memory leak
1269	//double saveY=new double[numberTotal];*
1270	//double saveS=new double[numberTotal];*
1271	//double saveY2=new double[numberTotal];*
1272	//double saveS2=new double[numberTotal];*
1273	double saveY[`100`];
1274	double saveS[`100`];
1275	double saveY2[`100`];
1276	double saveS2[`100`];
1277	double zz[`10000`];
1278	#endif
1279	double * dArray2 = dArray + numberTotal;
1280	// big big loop
1281	while (!allFinished) {
1282	double * work = longArray->denseVector();
1283	int * which = longArray->getIndices();
1284	allFinished = true;
1285	// CONJUGATE 0 - never, 1 as pivotMode, 2 as localPivotMode, 3 always
1286	//#define SMALLTHETA1 1.0e-25
1287	//#define SMALLTHETA2 1.0e-25
1288	#define SMALLTHETA1 1.0e-10
1289	#define SMALLTHETA2 1.0e-10
1290	#define CONJUGATE 2
1291	#if CONJUGATE == 0
1292	int conjugate = `0`;
1293	#elif CONJUGATE == 1
1294	int conjugate = (pivotMode < `10`) ? MINTYPE : `0`;
1295	#elif CONJUGATE == 2
1296	int conjugate = MINTYPE;
1297	#else
1298	int conjugate = MINTYPE;
1299	#endif
1300
1301	//if (pivotMode==1)
1302	//localPivotMode=11;
1303	#if CLP_DEBUG > 1
1304	longArray->checkClear();
1305	#endif
1306	int numberNonBasic = `0`;
1307	int startLocalMode = -`1`;
1308	while (!finished) {
1309	double simpleObjective = COIN_DBL_MAX;
1310	returnCode = `1`;
1311	int iSequence;
1312	objective_->reducedGradient(this, dj_, false);
1313	// get direction vector in longArray
1314	longArray->clear();
1315	// take out comment to try slightly different idea
1316	if (nPasses + nTotalPasses > `3000` \|\| nBigPasses > `100`) {
1317	if (factorization_->pivots())
1318	returnCode = `3`;
1319	break;
1320	}
1321	if (!nPasses) {
1322	numberNonBasic = `0`;
1323	nBigPasses++;
1324	}
1325	// just do superbasic if in cleanup mode
1326	int local = localPivotMode;
1327	if (!local && pivotMode >= `10` && nBigPasses < `10`) {
1328	local = `100`;
1329	} else if (local == -`1` \|\| nBigPasses >= `10`) {
1330	local = `0`;
1331	localPivotMode = `0`;
1332	}
1333	if (justOne) {
1334	local = `2`;
1335	//local=100;
1336	justOne = false;
1337	}
1338	if (!nPasses)
1339	startLocalMode = local;
1340	directionVector(longArray, spare, rowArray, local,
1341	normFlagged, normUnflagged, numberNonBasic);
1342	{
1343	// check null vector
1344	double * rhs = spare->denseVector();
1345	#if CLP_DEBUG > 1
1346	spare->checkClear();
1347	#endif
1348	int iRow;
1349	multiplyAdd(solution_ + numberColumns_, numberRows_, -`1.0`, rhs, `0.0`);
1350	matrix_->times(`1.0`, solution_, rhs, rowScale_, columnScale_);
1351	double largest = `0.0`;
1352	int iLargest = -`1`;
1353	for (iRow = `0`; iRow < numberRows_; iRow++) {
1354	double value = fabs(rhs[iRow]);
1355	rhs[iRow] = `0.0`;
1356	if (value > largest) {
1357	largest = value;
1358	iLargest = iRow;
1359	}
1360	}
1361	#if CLP_DEBUG > 0
1362	if ((handler_->logLevel() & `32`) && largest > `1.0e-8`)
1363	printf("largest rhs error %g on row %d\n", largest, iLargest);
1364	#endif
1365	if (solutionError < `0.0`) {
1366	solutionError = largest;
1367	} else if (largest > CoinMax(`1.0e-8`, `1.0e2` * solutionError) &&
1368	factorization_->pivots()) {
1369	longArray->clear();
1370	pivotRow_ = -`1`;
1371	theta_ = `0.0`;
1372	return `3`;
1373	}
1374	}
1375	if (sequenceIn_ >= `0`)
1376	lastSequenceIn = sequenceIn_;
1377	#if MINTYPE!=2
1378	double djNormSave = djNorm;
1379	#endif
1380	djNorm = `0.0`;
1381	int iIndex;
1382	for (iIndex = `0`; iIndex < numberNonBasic; iIndex++) {
1383	int iSequence = which[iIndex];
1384	double alpha = work[iSequence];
1385	djNorm += alpha * alpha;
1386	}
1387	// go to conjugate gradient if necessary
1388	if (numberNonBasic && localPivotMode >= `10` && (nPasses > `4` \|\| sequenceIn_ < `0`)) {
1389	localPivotMode = `0`;
1390	nTotalPasses += nPasses;
1391	nPasses = `0`;
1392	}
1393	#if CONJUGATE == 2
1394	conjugate = (localPivotMode < `10`) ? MINTYPE : `0`;
1395	#endif
1396	//printf("bigP %d pass %d nBasic %d norm %g normI %g normF %g\n",
1397	// nBigPasses,nPasses,numberNonBasic,normUnflagged,normFlagged);
1398	if (!nPasses) {
1399	//printf("numberNon %d\n",numberNonBasic);
1400	#if MINTYPE==2
1401	assert (numberNonBasic < `100`);
1402	memset(zz, `0`, numberNonBasic * numberNonBasic * sizeof(double));
1403	int put = `0`;
1404	for (int iVariable = `0`; iVariable < numberNonBasic; iVariable++) {
1405	zz[put] = `1.0`;
1406	put += numberNonBasic + `1`;
1407	}
1408	#endif
1409	djNorm0 = CoinMax(djNorm, `1.0e-20`);
1410	CoinMemcpyN(work, numberTotal, dArray);
1411	CoinMemcpyN(work, numberTotal, dArray2);
1412	if (sequenceIn_ >= `0` && numberNonBasic == `1`) {
1413	// see if simple move
1414	double objTheta2 = objective_->stepLength(this, solution_, work, `1.0e30`,
1415	currentObj, predictedObj, thetaObj);
1416	rowArray->clear();
1417
1418	// update the incoming column
1419	unpackPacked(rowArray);
1420	factorization_->updateColumnFT(spare, rowArray);
1421	theta_ = `0.0`;
1422	double * work2 = rowArray->denseVector();
1423	int number = rowArray->getNumElements();
1424	int * which2 = rowArray->getIndices();
1425	int iIndex;
1426	bool easyMove = false;
1427	double way;
1428	if (dj_[sequenceIn_] > `0.0`)
1429	way = -`1.0`;
1430	else
1431	way = `1.0`;
1432	double largest = COIN_DBL_MAX;
1433	int kPivot = -`1`;
1434	for (iIndex = `0`; iIndex < number; iIndex++) {
1435	int iRow = which2[iIndex];
1436	double alpha = way * work2[iIndex];
1437	int iPivot = pivotVariable_[iRow];
1438	if (alpha < -`1.0e-5`) {
1439	double distance = upper_[iPivot] - solution_[iPivot];
1440	if (distance < -largest * alpha) {
1441	kPivot = iPivot;
1442	largest = CoinMax(`0.0`, -distance / alpha);
1443	}
1444	if (distance < -`1.0e-12` * alpha) {
1445	easyMove = true;
1446	break;
1447	}
1448	} else if (alpha > `1.0e-5`) {
1449	double distance = solution_[iPivot] - lower_[iPivot];
1450	if (distance < largest * alpha) {
1451	kPivot = iPivot;
1452	largest = CoinMax(`0.0`, distance / alpha);
1453	}
1454	if (distance < `1.0e-12` * alpha) {
1455	easyMove = true;
1456	break;
1457	}
1458	}
1459	}
1460	// But largest has to match up with change
1461	assert (work[sequenceIn_]);
1462	largest /= fabs(work[sequenceIn_]);
1463	if (largest < objTheta2) {
1464	easyMove = true;
1465	} else if (!easyMove) {
1466	double objDrop = currentObj - predictedObj;
1467	double th = objective_->stepLength(this, solution_, work, largest,
1468	currentObj, predictedObj, simpleObjective);
1469	simpleObjective = CoinMax(simpleObjective, predictedObj);
1470	double easyDrop = currentObj - simpleObjective;
1471	if (easyDrop > `1.0e-8` && easyDrop > `0.5` * objDrop) {
1472	easyMove = true;
1473	#ifdef CLP_DEBUG
1474	if (handler_->logLevel() & `32`)
1475	printf("easy - obj drop %g, easy drop %g\n", objDrop, easyDrop);
1476	#endif
1477	if (easyDrop > objDrop) {
1478	// debug
1479	printf("****** th %g simple %g\n", th, simpleObjective);
1480	objective_->stepLength(this, solution_, work, `1.0e30`,
1481	currentObj, predictedObj, simpleObjective);
1482	objective_->stepLength(this, solution_, work, largest,
1483	currentObj, predictedObj, simpleObjective);
1484	}
1485	}
1486	}
1487	rowArray->clear();
1488	#ifdef CLP_DEBUG
1489	if (handler_->logLevel() & `32`)
1490	printf("largest %g piv %d\n", largest, kPivot);
1491	#endif
1492	if (easyMove) {
1493	valueIn_ = solution_[sequenceIn_];
1494	dualIn_ = dj_[sequenceIn_];
1495	lowerIn_ = lower_[sequenceIn_];
1496	upperIn_ = upper_[sequenceIn_];
1497	if (dualIn_ > `0.0`)
1498	directionIn_ = -`1`;
1499	else
1500	directionIn_ = `1`;
1501	longArray->clear();
1502	pivotRow_ = -`1`;
1503	theta_ = `0.0`;
1504	return `0`;
1505	}
1506	}
1507	} else {
1508	#if MINTYPE==1
1509	if (conjugate) {
1510	double djNorm2 = djNorm;
1511	#if 0
1512	if (numberNonBasic) {
1513	int iIndex;
1514	djNorm2 = `0.0`;
1515	for (iIndex = `0`; iIndex < numberNonBasic; iIndex++) {
1516	int iSequence = which[iIndex];
1517	double alpha = work[iSequence];
1518	//djNorm2 += alphaalpha;*
1519	double alpha2 = work[iSequence] - dArray2[iSequence];
1520	djNorm2 += alpha * alpha2;
1521	}
1522	//printf("a %.18g b %.18g\n",djNorm,djNorm2);
1523	}
1524	#endif
1525	djNorm = djNorm2;
1526	double beta = djNorm2 / djNormSave;
1527	// reset beta every so often
1528	//if (numberNonBasic&&nPasses>numberNonBasic&&(nPasses%(3numberNonBasic))==1)*
1529	//beta=0.0;
1530	//printf("current norm %g, old %g - beta %g\n",
1531	// djNorm,djNormSave,beta);
1532	for (iSequence = `0`; iSequence < numberTotal; iSequence++) {
1533	dArray[iSequence] = work[iSequence] + beta * dArray[iSequence];
1534	dArray2[iSequence] = work[iSequence];
1535	}
1536	} else {
1537	for (iSequence = `0`; iSequence < numberTotal; iSequence++)
1538	dArray[iSequence] = work[iSequence];
1539	}
1540	#else
1541	int number2 = numberNonBasic;
1542	if (number2) {
1543	// pack down into dArray
1544	int jLast = -`1`;
1545	for (iSequence = `0`; iSequence < numberNonBasic; iSequence++) {
1546	int j = which[iSequence];
1547	assert(j > jLast);
1548	jLast = j;
1549	double value = work[j];
1550	dArray[iSequence] = -value;
1551	}
1552	// see whether to restart
1553	// check signs - gradient
1554	double g1 = innerProduct(dArray, number2, dArray);
1555	double g2 = innerProduct(dArray, number2, saveY2);
1556	// Get differences
1557	for (iSequence = `0`; iSequence < numberNonBasic; iSequence++) {
1558	saveY2[iSequence] = dArray[iSequence] - saveY2[iSequence];
1559	saveS2[iSequence] = solution_[iSequence] - saveS2[iSequence];
1560	}
1561	double g3 = innerProduct(saveS2, number2, saveY2);
1562	printf("inner %g\n", g3);
1563	//assert(g3>0);
1564	double zzz[`10000`];
1565	int iVariable;
1566	g2 = `1.0e50`; // temp
1567	if (fabs(g2) >= `0.2` * fabs(g1)) {
1568	// restart
1569	double delta = innerProduct(saveY2, number2, saveS2) /
1570	innerProduct(saveY2, number2, saveY2);
1571	delta = `1.0`; //temp
1572	memset(zz, `0`, number2 * sizeof(double));
1573	int put = `0`;
1574	for (iVariable = `0`; iVariable < number2; iVariable++) {
1575	zz[put] = delta;
1576	put += number2 + `1`;
1577	}
1578	} else {
1579	}
1580	CoinMemcpyN(zz, number2 * number2, zzz);
1581	double ww[`100`];
1582	// get sk -Hkyk
1583	for (iVariable = `0`; iVariable < number2; iVariable++) {
1584	double value = `0.0`;
1585	for (int jVariable = `0`; jVariable < number2; jVariable++) {
1586	value += saveY2[jVariable] * zzz[iVariable+number2*jVariable];
1587	}
1588	ww[iVariable] = saveS2[iVariable] - value;
1589	}
1590	double ys = innerProduct(saveY2, number2, saveS2);
1591	double multiplier1 = `1.0` / ys;
1592	double multiplier2 = innerProduct(saveY2, number2, ww) / (ys * ys);
1593	#if 1
1594	// and second way
1595	// Hy
1596	double h[`100`];
1597	for (iVariable = `0`; iVariable < number2; iVariable++) {
1598	double value = `0.0`;
1599	for (int jVariable = `0`; jVariable < number2; jVariable++) {
1600	value += saveY2[jVariable] * zzz[iVariable+number2*jVariable];
1601	}
1602	h[iVariable] = value;
1603	}
1604	double hh[`10000`];
1605	double yhy1 = innerProduct(h, number2, saveY2) * multiplier1 + `1.0`;
1606	yhy1 *= multiplier1;
1607	for (iVariable = `0`; iVariable < number2; iVariable++) {
1608	for (int jVariable = `0`; jVariable < number2; jVariable++) {
1609	int put = iVariable + number2 * jVariable;
1610	double value = zzz[put];
1611	value += yhy1 * saveS2[iVariable] * saveS2[jVariable];
1612	hh[put] = value;
1613	}
1614	}
1615	for (iVariable = `0`; iVariable < number2; iVariable++) {
1616	for (int jVariable = `0`; jVariable < number2; jVariable++) {
1617	int put = iVariable + number2 * jVariable;
1618	double value = hh[put];
1619	value -= multiplier1 * (saveS2[iVariable] * h[jVariable]);
1620	value -= multiplier1 * (saveS2[jVariable] * h[iVariable]);
1621	hh[put] = value;
1622	}
1623	}
1624	#endif
1625	// now update H
1626	for (iVariable = `0`; iVariable < number2; iVariable++) {
1627	for (int jVariable = `0`; jVariable < number2; jVariable++) {
1628	int put = iVariable + number2 * jVariable;
1629	double value = zzz[put];
1630	value += multiplier1 * (ww[iVariable] * saveS2[jVariable]
1631	+ ww[jVariable] * saveS2[iVariable]);
1632	value -= multiplier2 * saveS2[iVariable] * saveS2[jVariable];
1633	zzz[put] = value;
1634	}
1635	}
1636	//memcpy(zzz,hh,sizesizeof(double));*
1637	// do search direction
1638	memset(dArray, `0`, numberTotal * sizeof(double));
1639	for (iVariable = `0`; iVariable < numberNonBasic; iVariable++) {
1640	double value = `0.0`;
1641	for (int jVariable = `0`; jVariable < number2; jVariable++) {
1642	int k = which[jVariable];
1643	value += work[k] * zzz[iVariable+number2*jVariable];
1644	}
1645	int i = which[iVariable];
1646	dArray[i] = value;
1647	}
1648	// Now fill out dArray
1649	{
1650	int j;
1651	// Now do basic ones
1652	int iRow;
1653	CoinIndexedVector * spare1 = spare;
1654	CoinIndexedVector * spare2 = rowArray;
1655	#if CLP_DEBUG > 1
1656	spare1->checkClear();
1657	spare2->checkClear();
1658	#endif
1659	double * array2 = spare1->denseVector();
1660	int * index2 = spare1->getIndices();
1661	int number2 = `0`;
1662	times(-`1.0`, dArray, array2);
1663	dArray = dArray + numberColumns_;
1664	for (iRow = `0`; iRow < numberRows_; iRow++) {
1665	double value = array2[iRow] + dArray[iRow];
1666	if (value) {
1667	array2[iRow] = value;
1668	index2[number2++] = iRow;
1669	} else {
1670	array2[iRow] = `0.0`;
1671	}
1672	}
1673	dArray -= numberColumns_;
1674	spare1->setNumElements(number2);
1675	// Ftran
1676	factorization_->updateColumn(spare2, spare1);
1677	number2 = spare1->getNumElements();
1678	for (j = `0`; j < number2; j++) {
1679	int iSequence = index2[j];
1680	double value = array2[iSequence];
1681	array2[iSequence] = `0.0`;
1682	if (value) {
1683	int iPivot = pivotVariable_[iSequence];
1684	double oldValue = dArray[iPivot];
1685	dArray[iPivot] = value + oldValue;
1686	}
1687	}
1688	spare1->setNumElements(`0`);
1689	}
1690	//assert (innerProductIndexed(dArray,number2,work,which)>0);
1691	//printf ("innerP1 %g\n",innerProduct(dArray,numberTotal,work));
1692	printf ("innerP1 %g innerP2 %g\n", innerProduct(dArray, numberTotal, work),
1693	innerProductIndexed(dArray, numberNonBasic, work, which));
1694	assert (innerProduct(dArray, numberTotal, work) > `0`);
1695	#if 1
1696	{
1697	// check null vector
1698	int iRow;
1699	double qq[`10`];
1700	memset(qq, `0`, numberRows_ * sizeof(double));
1701	multiplyAdd(dArray + numberColumns_, numberRows_, -`1.0`, qq, `0.0`);
1702	matrix_->times(`1.0`, dArray, qq, rowScale_, columnScale_);
1703	double largest = `0.0`;
1704	int iLargest = -`1`;
1705	for (iRow = `0`; iRow < numberRows_; iRow++) {
1706	double value = fabs(qq[iRow]);
1707	if (value > largest) {
1708	largest = value;
1709	iLargest = iRow;
1710	}
1711	}
1712	printf("largest null error %g on row %d\n", largest, iLargest);
1713	for (iSequence = `0`; iSequence < numberTotal; iSequence++) {
1714	if (getStatus(iSequence) == basic)
1715	assert (fabs(dj_[iSequence]) < `1.0e-3`);
1716	}
1717	}
1718	#endif
1719	CoinMemcpyN(saveY2, numberNonBasic, saveY);
1720	CoinMemcpyN(saveS2, numberNonBasic, saveS);
1721	}
1722	#endif
1723	}
1724	#if MINTYPE==2
1725	for (iSequence = `0`; iSequence < numberNonBasic; iSequence++) {
1726	int j = which[iSequence];
1727	saveY2[iSequence] = -work[j];
1728	saveS2[iSequence] = solution_[j];
1729	}
1730	#endif
1731	if (djNorm < eps * djNorm0 \|\| (nPasses > `100` && djNorm < CoinMin(`1.0e-1` * djNorm0, `1.0e-12`))) {
1732	#ifdef CLP_DEBUG
1733	if (handler_->logLevel() & `32`)
1734	printf("dj norm reduced from %g to %g\n", djNorm0, djNorm);
1735	#endif
1736	if (pivotMode < `10` \|\| !numberNonBasic) {
1737	finished = true;
1738	} else {
1739	localPivotMode = pivotMode;
1740	nTotalPasses += nPasses;
1741	nPasses = `0`;
1742	continue;
1743	}
1744	}
1745	//if (nPasses==100\|\|nPasses==50)
1746	//printf("pass %d dj norm reduced from %g to %g - flagged norm %g\n",nPasses,djNorm0,djNorm,
1747	// normFlagged);
1748	if (nPasses > `100` && djNorm < `1.0e-2` * normFlagged && !startLocalMode) {
1749	#ifdef CLP_DEBUG
1750	if (handler_->logLevel() & `32`)
1751	printf("dj norm reduced from %g to %g - flagged norm %g - unflagging\n", djNorm0, djNorm,
1752	normFlagged);
1753	#endif
1754	unflag();
1755	localPivotMode = `0`;
1756	nTotalPasses += nPasses;
1757	nPasses = `0`;
1758	continue;
1759	}
1760	if (djNorm > `0.99` * djNorm0 && nPasses > `1500`) {
1761	finished = true;
1762	#ifdef CLP_DEBUG
1763	if (handler_->logLevel() & `32`)
1764	printf("dj norm NOT reduced from %g to %g\n", djNorm0, djNorm);
1765	#endif
1766	djNorm = `1.2345e-20`;
1767	}
1768	number = longArray->getNumElements();
1769	if (!numberNonBasic) {
1770	// looks optimal
1771	// check if any dj look plausible
1772	int nSuper = `0`;
1773	int nFlagged = `0`;
1774	int nNormal = `0`;
1775	for (int iSequence = `0`; iSequence < numberColumns_ + numberRows_; iSequence++) {
1776	if (flagged(iSequence)) {
1777	switch(getStatus(iSequence)) {
1778
1779	case basic:
1780	case ClpSimplex::isFixed:
1781	break;
1782	case atUpperBound:
1783	if (dj_[iSequence] > dualTolerance_)
1784	nFlagged++;
1785	break;
1786	case atLowerBound:
1787	if (dj_[iSequence] < -dualTolerance_)
1788	nFlagged++;
1789	break;
1790	case isFree:
1791	case superBasic:
1792	if (fabs(dj_[iSequence]) > dualTolerance_)
1793	nFlagged++;
1794	break;
1795	}
1796	continue;
1797	}
1798	switch(getStatus(iSequence)) {
1799
1800	case basic:
1801	case ClpSimplex::isFixed:
1802	break;
1803	case atUpperBound:
1804	if (dj_[iSequence] > dualTolerance_)
1805	nNormal++;
1806	break;
1807	case atLowerBound:
1808	if (dj_[iSequence] < -dualTolerance_)
1809	nNormal++;
1810	break;
1811	case isFree:
1812	case superBasic:
1813	if (fabs(dj_[iSequence]) > dualTolerance_)
1814	nSuper++;
1815	break;
1816	}
1817	}
1818	#ifdef CLP_DEBUG
1819	if (handler_->logLevel() & `32`)
1820	printf("%d super, %d normal, %d flagged\n",
1821	nSuper, nNormal, nFlagged);
1822	#endif
1823
1824	int nFlagged2 = `1`;
1825	if (lastSequenceIn < `0` && !nNormal && !nSuper) {
1826	nFlagged2 = unflag();
1827	if (pivotMode >= `10`) {
1828	pivotMode--;
1829	#ifdef CLP_DEBUG
1830	if (handler_->logLevel() & `32`)
1831	printf("pivot mode now %d\n", pivotMode);
1832	#endif
1833	if (pivotMode == `9`)
1834	pivotMode = `0`; // switch off fast attempt
1835	}
1836	} else {
1837	lastSequenceIn = -`1`;
1838	}
1839	if (!nFlagged2 \|\| (normFlagged + normUnflagged < `1.0e-8`)) {
1840	primalColumnPivot_->setLooksOptimal(true);
1841	return `0`;
1842	} else {
1843	localPivotMode = -`1`;
1844	nTotalPasses += nPasses;
1845	nPasses = `0`;
1846	finished = false;
1847	continue;
1848	}
1849	}
1850	bestSequence = -`1`;
1851	bestBasicSequence = -`1`;
1852
1853	// temp
1854	nPasses++;
1855	if (nPasses > `2000`)
1856	finished = true;
1857	double theta = `1.0e30`;
1858	basicTheta = `1.0e30`;
1859	theta_ = `1.0e30`;
1860	double basicTolerance = `1.0e-4` * primalTolerance_;
1861	for (iSequence = `0`; iSequence < numberTotal; iSequence++) {
1862	//if (flagged(iSequence)
1863	// continue;
1864	double alpha = dArray[iSequence];
1865	Status thisStatus = getStatus(iSequence);
1866	double oldValue = solution_[iSequence];
1867	if (thisStatus != basic) {
1868	if (fabs(alpha) >= acceptablePivot) {
1869	if (alpha < `0.0`) {
1870	// variable going towards lower bound
1871	double bound = lower_[iSequence];
1872	oldValue -= bound;
1873	if (oldValue + theta * alpha < `0.0`) {
1874	bestSequence = iSequence;
1875	theta = CoinMax(`0.0`, oldValue / (-alpha));
1876	}
1877	} else {
1878	// variable going towards upper bound
1879	double bound = upper_[iSequence];
1880	oldValue = bound - oldValue;
1881	if (oldValue - theta * alpha < `0.0`) {
1882	bestSequence = iSequence;
1883	theta = CoinMax(`0.0`, oldValue / alpha);
1884	}
1885	}
1886	}
1887	} else {
1888	if (fabs(alpha) >= acceptableBasic) {
1889	if (alpha < `0.0`) {
1890	// variable going towards lower bound
1891	double bound = lower_[iSequence];
1892	oldValue -= bound;
1893	if (oldValue + basicTheta * alpha < -basicTolerance) {
1894	bestBasicSequence = iSequence;
1895	basicTheta = CoinMax(`0.0`, (oldValue + basicTolerance) / (-alpha));
1896	}
1897	} else {
1898	// variable going towards upper bound
1899	double bound = upper_[iSequence];
1900	oldValue = bound - oldValue;
1901	if (oldValue - basicTheta * alpha < -basicTolerance) {
1902	bestBasicSequence = iSequence;
1903	basicTheta = CoinMax(`0.0`, (oldValue + basicTolerance) / alpha);
1904	}
1905	}
1906	}
1907	}
1908	}
1909	theta_ = CoinMin(theta, basicTheta);
1910	// Now find minimum of function
1911	double objTheta2 = objective_->stepLength(this, solution_, dArray, CoinMin(theta, basicTheta),
1912	currentObj, predictedObj, thetaObj);
1913	#ifdef CLP_DEBUG
1914	if (handler_->logLevel() & `32`)
1915	printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj);
1916	#endif
1917	#if MINTYPE==1
1918	if (conjugate) {
1919	double offset;
1920	const double * gradient = objective_->gradient(this,
1921	dArray, offset,
1922	true, `0`);
1923	double product = `0.0`;
1924	for (iSequence = `0`; iSequence < numberColumns_; iSequence++) {
1925	double alpha = dArray[iSequence];
1926	double value = alpha * gradient[iSequence];
1927	product += value;
1928	}
1929	//#define INCLUDESLACK
1930	#ifdef INCLUDESLACK
1931	for (; iSequence < numberColumns_ + numberRows_; iSequence++) {
1932	double alpha = dArray[iSequence];
1933	double value = alpha * cost_[iSequence];
1934	product += value;
1935	}
1936	#endif
1937	if (product > `0.0`)
1938	objTheta = djNorm / product;
1939	else
1940	objTheta = COIN_DBL_MAX;
1941	#ifndef NDEBUG
1942	if (product < -`1.0e-8` && handler_->logLevel() > `1`)
1943	printf("bad product %g\n", product);
1944	#endif
1945	product = CoinMax(product, `0.0`);
1946	} else {
1947	objTheta = objTheta2;
1948	}
1949	#else
1950	objTheta = objTheta2;
1951	#endif
1952	// if very small difference then take pivot (depends on djNorm?)
1953	// distinguish between basic and non-basic
1954	bool chooseObjTheta = objTheta < theta_;
1955	if (chooseObjTheta) {
1956	if (thetaObj < currentObj - `1.0e-12` && objTheta + `1.0e-10` > theta_)
1957	chooseObjTheta = false;
1958	//if (thetaObj<currentObj+1.0e-12&&objTheta+1.0e-5>theta_)
1959	//chooseObjTheta=false;
1960	}
1961	//if (objTheta+SMALLTHETA1<theta_\|\|(thetaObj>currentObj+difference&&objTheta<theta_)) {
1962	if (chooseObjTheta) {
1963	theta_ = objTheta;
1964	} else {
1965	objTheta = CoinMax(objTheta, `1.00000001` * theta_ + `1.0e-12`);
1966	//if (theta+1.0e-13>basicTheta) {
1967	//theta = CoinMax(theta,1.00000001basicTheta);*
1968	//theta_ = basicTheta;
1969	//}
1970	}
1971	// always out if one variable in and zero move
1972	if (theta_ == basicTheta \|\| (sequenceIn_ >= `0` && theta_ < `1.0e-10`))
1973	finished = true; // come out
1974	#ifdef CLP_DEBUG
1975	if (handler_->logLevel() & `32`) {
1976	printf("%d pass,", nPasses);
1977	if (sequenceIn_ >= `0`)
1978	printf (" Sin %d,", sequenceIn_);
1979	if (basicTheta == theta_)
1980	printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta);
1981	else
1982	printf(" basicTheta %g", basicTheta);
1983	if (theta == theta_)
1984	printf(" X(%d) non-basicTheta %g", bestSequence, theta);
1985	else
1986	printf(" non-basicTheta %g", theta);
1987	printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta);
1988	printf(" djNorm %g\n", djNorm);
1989	}
1990	#endif
1991	if (handler_->logLevel() > `3` && objTheta != theta_) {
1992	printf("%d pass obj %g,", nPasses, currentObj);
1993	if (sequenceIn_ >= `0`)
1994	printf (" Sin %d,", sequenceIn_);
1995	if (basicTheta == theta_)
1996	printf(" X(%d) basicTheta %g", bestBasicSequence, basicTheta);
1997	else
1998	printf(" basicTheta %g", basicTheta);
1999	if (theta == theta_)
2000	printf(" X(%d) non-basicTheta %g", bestSequence, theta);
2001	else
2002	printf(" non-basicTheta %g", theta);
2003	printf(" %s objTheta %g", objTheta == theta_ ? "X" : "", objTheta);
2004	printf(" djNorm %g\n", djNorm);
2005	}
2006	if (objTheta != theta_) {
2007	//printf("hit boundary after %d passes\n",nPasses);
2008	nTotalPasses += nPasses;
2009	nPasses = `0`; // start again
2010	}
2011	if (localPivotMode < `10` \|\| lastSequenceIn == bestSequence) {
2012	if (theta_ == theta && theta_ < basicTheta && theta_ < `1.0e-5`)
2013	setFlagged(bestSequence); // out of active set
2014	}
2015	// Update solution
2016	for (iSequence = `0`; iSequence < numberTotal; iSequence++) {
2017	//for (iIndex=0;iIndex<number;iIndex++) {
2018
2019	//int iSequence = which[iIndex];
2020	double alpha = dArray[iSequence];
2021	if (alpha) {
2022	double value = solution_[iSequence] + theta_ * alpha;
2023	solution_[iSequence] = value;
2024	switch(getStatus(iSequence)) {
2025
2026	case basic:
2027	case isFixed:
2028	case isFree:
2029	case atUpperBound:
2030	case atLowerBound:
2031	nonLinearCost_->setOne(iSequence, value);
2032	break;
2033	case superBasic:
2034	// To get correct action
2035	setStatus(iSequence, isFixed);
2036	nonLinearCost_->setOne(iSequence, value);
2037	//assert (getStatus(iSequence)!=isFixed);
2038	break;
2039	}
2040	}
2041	}
2042	if (objTheta < SMALLTHETA2 && objTheta == theta_) {
2043	if (sequenceIn_ >= `0` && basicTheta < `1.0e-9`) {
2044	// try just one
2045	localPivotMode = `0`;
2046	sequenceIn_ = -`1`;
2047	nTotalPasses += nPasses;
2048	nPasses = `0`;
2049	//finished=true;
2050	//objTheta=0.0;
2051	//theta_=0.0;
2052	} else if (sequenceIn_ < `0` && nTotalPasses > `10`) {
2053	if (objTheta < `1.0e-10`) {
2054	finished = true;
2055	//printf("zero move\n");
2056	break;
2057	}
2058	}
2059	}
2060	#ifdef CLP_DEBUG
2061	if (handler_->logLevel() & `32`) {
2062	objective_->stepLength(this, solution_, work, `0.0`,
2063	currentObj, predictedObj, thetaObj);
2064	printf("current obj %g after update - simple was %g\n", currentObj, simpleObjective);
2065	}
2066	#endif
2067	if (sequenceIn_ >= `0` && !finished && objTheta > `1.0e-4`) {
2068	// we made some progress - back to normal
2069	if (localPivotMode < `10`) {
2070	localPivotMode = `0`;
2071	sequenceIn_ = -`1`;
2072	nTotalPasses += nPasses;
2073	nPasses = `0`;
2074	}
2075	#ifdef CLP_DEBUG
2076	if (handler_->logLevel() & `32`)
2077	printf("some progress?\n");
2078	#endif
2079	}
2080	#if CLP_DEBUG > 1
2081	longArray->checkClean();
2082	#endif
2083	}
2084	#ifdef CLP_DEBUG
2085	if (handler_->logLevel() & `32`)
2086	printf("out of loop after %d (%d) passes\n", nPasses, nTotalPasses);
2087	#endif
2088	if (nTotalPasses >= `1000` \|\| (nTotalPasses > `10` && sequenceIn_ < `0` && theta_ < `1.0e-10`))
2089	returnCode = `2`;
2090	bool ordinaryDj = false;
2091	//if(sequenceIn_>=0&&numberNonBasic==1&&theta_<1.0e-7&&theta_==basicTheta)
2092	//printf("could try ordinary iteration %g\n",theta_);
2093	if(sequenceIn_ >= `0` && numberNonBasic == `1` && theta_ < `1.0e-15`) {
2094	//printf("trying ordinary iteration\n");
2095	ordinaryDj = true;
2096	}
2097	if (!basicTheta && !ordinaryDj) {
2098	//returnCode=2;
2099	//objTheta=-1.0; // so we fall through
2100	}
2101	assert (theta_ < `1.0e30`); // for now
2102	// See if we need to pivot
2103	if (theta_ == basicTheta \|\| ordinaryDj) {
2104	if (!ordinaryDj) {
2105	// find an incoming column which will force pivot
2106	int iRow;
2107	pivotRow_ = -`1`;
2108	for (iRow = `0`; iRow < numberRows_; iRow++) {
2109	if (pivotVariable_[iRow] == bestBasicSequence) {
2110	pivotRow_ = iRow;
2111	break;
2112	}
2113	}
2114	assert (pivotRow_ >= `0`);
2115	// Get good size for pivot
2116	double acceptablePivot = `1.0e-7`;
2117	if (factorization_->pivots() > `10`)
2118	acceptablePivot = `1.0e-5`; // if we have iterated be more strict
2119	else if (factorization_->pivots() > `5`)
2120	acceptablePivot = `1.0e-6`; // if we have iterated be slightly more strict
2121	// should be dArray but seems better this way!
2122	double direction = work[bestBasicSequence] > `0.0` ? -`1.0` : `1.0`;
2123	// create as packed
2124	rowArray->createPacked(`1`, &pivotRow_, &direction);
2125	factorization_->updateColumnTranspose(spare, rowArray);
2126	// put row of tableau in rowArray and columnArray
2127	matrix_->transposeTimes(this, -`1.0`,
2128	rowArray, spare, columnArray);
2129	// choose one futhest away from bound which has reasonable pivot
2130	// If increasing we want negative alpha
2131
2132	double * work2;
2133	int iSection;
2134
2135	sequenceIn_ = -`1`;
2136	double bestValue = -`1.0`;
2137	double bestDirection = `0.0`;
2138	// First pass we take correct direction and large pivots
2139	// then correct direction
2140	// then any
2141	double check[] = {`1.0e-8`, -`1.0e-12`, -`1.0e30`};
2142	double mult[] = {`100.0`, `1.0`, `1.0`};
2143	for (int iPass = `0`; iPass < `3`; iPass++) {
2144	//if (!bestValue&&iPass==2)
2145	//bestValue=-1.0;
2146	double acceptable = acceptablePivot * mult[iPass];
2147	double checkValue = check[iPass];
2148	for (iSection = `0`; iSection < `2`; iSection++) {
2149
2150	int addSequence;
2151
2152	if (!iSection) {
2153	work2 = rowArray->denseVector();
2154	number = rowArray->getNumElements();
2155	which = rowArray->getIndices();
2156	addSequence = numberColumns_;
2157	} else {
2158	work2 = columnArray->denseVector();
2159	number = columnArray->getNumElements();
2160	which = columnArray->getIndices();
2161	addSequence = `0`;
2162	}
2163	int i;
2164
2165	for (i = `0`; i < number; i++) {
2166	int iSequence = which[i] + addSequence;
2167	if (flagged(iSequence))
2168	continue;
2169	//double distance = CoinMin(solution_[iSequence]-lower_[iSequence],
2170	// upper_[iSequence]-solution_[iSequence]);
2171	double alpha = work2[i];
2172	// should be dArray but seems better this way!
2173	double change = work[iSequence];
2174	Status thisStatus = getStatus(iSequence);
2175	double direction = `0`;
2176	switch(thisStatus) {
2177
2178	case basic:
2179	case ClpSimplex::isFixed:
2180	break;
2181	case isFree:
2182	case superBasic:
2183	if (alpha < -acceptable && change > checkValue)
2184	direction = `1.0`;
2185	else if (alpha > acceptable && change < -checkValue)
2186	direction = -`1.0`;
2187	break;
2188	case atUpperBound:
2189	if (alpha > acceptable && change < -checkValue)
2190	direction = -`1.0`;
2191	else if (iPass == `2` && alpha < -acceptable && change < -checkValue)
2192	direction = `1.0`;
2193	break;
2194	case atLowerBound:
2195	if (alpha < -acceptable && change > checkValue)
2196	direction = `1.0`;
2197	else if (iPass == `2` && alpha > acceptable && change > checkValue)
2198	direction = -`1.0`;
2199	break;
2200	}
2201	if (direction) {
2202	if (sequenceIn_ != lastSequenceIn \|\| localPivotMode < `10`) {
2203	if (CoinMin(solution_[iSequence] - lower_[iSequence],
2204	upper_[iSequence] - solution_[iSequence]) > bestValue) {
2205	bestValue = CoinMin(solution_[iSequence] - lower_[iSequence],
2206	upper_[iSequence] - solution_[iSequence]);
2207	sequenceIn_ = iSequence;
2208	bestDirection = direction;
2209	}
2210	} else {
2211	// choose
2212	bestValue = COIN_DBL_MAX;
2213	sequenceIn_ = iSequence;
2214	bestDirection = direction;
2215	}
2216	}
2217	}
2218	}
2219	if (sequenceIn_ >= `0` && bestValue > `0.0`)
2220	break;
2221	}
2222	if (sequenceIn_ >= `0`) {
2223	valueIn_ = solution_[sequenceIn_];
2224	dualIn_ = dj_[sequenceIn_];
2225	if (bestDirection < `0.0`) {
2226	// we want positive dj
2227	if (dualIn_ <= `0.0`) {
2228	//printf("bad dj - xx %g\n",dualIn_);
2229	dualIn_ = `1.0`;
2230	}
2231	} else {
2232	// we want negative dj
2233	if (dualIn_ >= `0.0`) {
2234	//printf("bad dj - xx %g\n",dualIn_);
2235	dualIn_ = -`1.0`;
2236	}
2237	}
2238	lowerIn_ = lower_[sequenceIn_];
2239	upperIn_ = upper_[sequenceIn_];
2240	if (dualIn_ > `0.0`)
2241	directionIn_ = -`1`;
2242	else
2243	directionIn_ = `1`;
2244	} else {
2245	//ordinaryDj=true;
2246	#ifdef CLP_DEBUG
2247	if (handler_->logLevel() & `32`) {
2248	printf("no easy pivot - norm %g mode %d\n", djNorm, localPivotMode);
2249	if (rowArray->getNumElements() + columnArray->getNumElements() < `12`) {
2250	for (iSection = `0`; iSection < `2`; iSection++) {
2251
2252	int addSequence;
2253
2254	if (!iSection) {
2255	work2 = rowArray->denseVector();
2256	number = rowArray->getNumElements();
2257	which = rowArray->getIndices();
2258	addSequence = numberColumns_;
2259	} else {
2260	work2 = columnArray->denseVector();
2261	number = columnArray->getNumElements();
2262	which = columnArray->getIndices();
2263	addSequence = `0`;
2264	}
2265	int i;
2266	char section[] = {`'R'`, `'C'`};
2267	for (i = `0`; i < number; i++) {
2268	int iSequence = which[i] + addSequence;
2269	if (flagged(iSequence)) {
2270	printf("%c%d flagged\n", section[iSection], which[i]);
2271	continue;
2272	} else {
2273	printf("%c%d status %d sol %g %g %g alpha %g change %g\n",
2274	section[iSection], which[i], status_[iSequence],
2275	lower_[iSequence], solution_[iSequence], upper_[iSequence],
2276	work2[i], work[iSequence]);
2277	}
2278	}
2279	}
2280	}
2281	}
2282	#endif
2283	assert (sequenceIn_ < `0`);
2284	justOne = true;
2285	allFinished = false; // Round again
2286	finished = false;
2287	nTotalPasses += nPasses;
2288	nPasses = `0`;
2289	if (djNorm < `0.9` * djNorm0 && djNorm < `1.0e-3` && !localPivotMode) {
2290	#ifdef CLP_DEBUG
2291	if (handler_->logLevel() & `32`)
2292	printf("no pivot - mode %d norms %g %g - unflagging\n",
2293	localPivotMode, djNorm0, djNorm);
2294	#endif
2295	unflag(); //unflagging
2296	returnCode = `1`;
2297	} else {
2298	returnCode = `2`; // do single incoming
2299	returnCode = `1`;
2300	}
2301	}
2302	}
2303	rowArray->clear();
2304	columnArray->clear();
2305	longArray->clear();
2306	if (ordinaryDj) {
2307	valueIn_ = solution_[sequenceIn_];
2308	dualIn_ = dj_[sequenceIn_];
2309	lowerIn_ = lower_[sequenceIn_];
2310	upperIn_ = upper_[sequenceIn_];
2311	if (dualIn_ > `0.0`)
2312	directionIn_ = -`1`;
2313	else
2314	directionIn_ = `1`;
2315	}
2316	if (returnCode == `1`)
2317	returnCode = `0`;
2318	} else {
2319	// round again
2320	longArray->clear();
2321	if (djNorm < `1.0e-3` && !localPivotMode) {
2322	if (djNorm == `1.2345e-20` && djNorm0 > `1.0e-4`) {
2323	#ifdef CLP_DEBUG
2324	if (handler_->logLevel() & `32`)
2325	printf("slow convergence djNorm0 %g, %d passes, mode %d, result %d\n", djNorm0, nPasses,
2326	localPivotMode, returnCode);
2327	#endif
2328	//if (!localPivotMode)
2329	//returnCode=2; // force singleton
2330	} else {
2331	#ifdef CLP_DEBUG
2332	if (handler_->logLevel() & `32`)
2333	printf("unflagging as djNorm %g %g, %d passes\n", djNorm, djNorm0, nPasses);
2334	#endif
2335	if (pivotMode >= `10`) {
2336	pivotMode--;
2337	#ifdef CLP_DEBUG
2338	if (handler_->logLevel() & `32`)
2339	printf("pivot mode now %d\n", pivotMode);
2340	#endif
2341	if (pivotMode == `9`)
2342	pivotMode = `0`; // switch off fast attempt
2343	}
2344	bool unflagged = unflag() != `0`;
2345	if (!unflagged && djNorm < `1.0e-10`) {
2346	// ? declare victory
2347	sequenceIn_ = -`1`;
2348	returnCode = `0`;
2349	}
2350	}
2351	}
2352	}
2353	}
2354	if (djNorm0 < `1.0e-12` * normFlagged) {
2355	#ifdef CLP_DEBUG
2356	if (handler_->logLevel() & `32`)
2357	printf("unflagging as djNorm %g %g and flagged norm %g\n", djNorm, djNorm0, normFlagged);
2358	#endif
2359	unflag();
2360	}
2361	if (saveObj - currentObj < `1.0e-5` && nTotalPasses > `2000`) {
2362	normUnflagged = `0.0`;
2363	double dualTolerance3 = CoinMin(`1.0e-2`, `1.0e3` * dualTolerance_);
2364	for (int iSequence = `0`; iSequence < numberColumns_ + numberRows_; iSequence++) {
2365	switch(getStatus(iSequence)) {
2366
2367	case basic:
2368	case ClpSimplex::isFixed:
2369	break;
2370	case atUpperBound:
2371	if (dj_[iSequence] > dualTolerance3)
2372	normFlagged += dj_[iSequence] * dj_[iSequence];
2373	break;
2374	case atLowerBound:
2375	if (dj_[iSequence] < -dualTolerance3)
2376	normFlagged += dj_[iSequence] * dj_[iSequence];
2377	break;
2378	case isFree:
2379	case superBasic:
2380	if (fabs(dj_[iSequence]) > dualTolerance3)
2381	normFlagged += dj_[iSequence] * dj_[iSequence];
2382	break;
2383	}
2384	}
2385	if (handler_->logLevel() > `2`)
2386	printf("possible optimal %d %d %g %g\n",
2387	nBigPasses, nTotalPasses, saveObj - currentObj, normFlagged);
2388	if (normFlagged < `1.0e-5`) {
2389	unflag();
2390	primalColumnPivot_->setLooksOptimal(true);
2391	returnCode = `0`;
2392	}
2393	}
2394	return returnCode;
2395	}
2396	/ Do last half of an iteration.*
2397	Return codes
2398	Reasons to come out normal mode
2399	-1 normal
2400	-2 factorize now - good iteration
2401	-3 slight inaccuracy - refactorize - iteration done
2402	-4 inaccuracy - refactorize - no iteration
2403	-5 something flagged - go round again
2404	+2 looks unbounded
2405	+3 max iterations (iteration done)
2406
2407	*/
2408	int
2409	ClpSimplexNonlinear::pivotNonlinearResult()
2410	{
2411
2412	int returnCode = -`1`;
2413
2414	rowArray_[`1`]->clear();
2415
2416	// we found a pivot column
2417	// update the incoming column
2418	unpackPacked(rowArray_[`1`]);
2419	factorization_->updateColumnFT(rowArray_[`2`], rowArray_[`1`]);
2420	theta_ = `0.0`;
2421	double * work = rowArray_[`1`]->denseVector();
2422	int number = rowArray_[`1`]->getNumElements();
2423	int * which = rowArray_[`1`]->getIndices();
2424	bool keepValue = false;
2425	double saveValue = `0.0`;
2426	if (pivotRow_ >= `0`) {
2427	sequenceOut_ = pivotVariable_[pivotRow_];
2428	valueOut_ = solution(sequenceOut_);
2429	keepValue = true;
2430	saveValue = valueOut_;
2431	lowerOut_ = lower_[sequenceOut_];
2432	upperOut_ = upper_[sequenceOut_];
2433	int iIndex;
2434	for (iIndex = `0`; iIndex < number; iIndex++) {
2435
2436	int iRow = which[iIndex];
2437	if (iRow == pivotRow_) {
2438	alpha_ = work[iIndex];
2439	break;
2440	}
2441	}
2442	} else {
2443	int iIndex;
2444	double smallest = COIN_DBL_MAX;
2445	for (iIndex = `0`; iIndex < number; iIndex++) {
2446	int iRow = which[iIndex];
2447	double alpha = work[iIndex];
2448	if (fabs(alpha) > `1.0e-6`) {
2449	int iPivot = pivotVariable_[iRow];
2450	double distance = CoinMin(upper_[iPivot] - solution_[iPivot],
2451	solution_[iPivot] - lower_[iPivot]);
2452	if (distance < smallest) {
2453	pivotRow_ = iRow;
2454	alpha_ = alpha;
2455	smallest = distance;
2456	}
2457	}
2458	}
2459	if (smallest > primalTolerance_) {
2460	smallest = COIN_DBL_MAX;
2461	for (iIndex = `0`; iIndex < number; iIndex++) {
2462	int iRow = which[iIndex];
2463	double alpha = work[iIndex];
2464	if (fabs(alpha) > `1.0e-6`) {
2465	double distance = randomNumberGenerator_.randomDouble();
2466	if (distance < smallest) {
2467	pivotRow_ = iRow;
2468	alpha_ = alpha;
2469	smallest = distance;
2470	}
2471	}
2472	}
2473	}
2474	assert (pivotRow_ >= `0`);
2475	sequenceOut_ = pivotVariable_[pivotRow_];
2476	valueOut_ = solution(sequenceOut_);
2477	lowerOut_ = lower_[sequenceOut_];
2478	upperOut_ = upper_[sequenceOut_];
2479	}
2480	double newValue = valueOut_ - theta_ * alpha_;
2481	bool isSuperBasic = false;
2482	if (valueOut_ >= upperOut_ - primalTolerance_) {
2483	directionOut_ = -`1`; // to upper bound
2484	upperOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
2485	upperOut_ = newValue;
2486	} else if (valueOut_ <= lowerOut_ + primalTolerance_) {
2487	directionOut_ = `1`; // to lower bound
2488	lowerOut_ = nonLinearCost_->nearest(sequenceOut_, newValue);
2489	} else {
2490	lowerOut_ = valueOut_;
2491	upperOut_ = valueOut_;
2492	isSuperBasic = true;
2493	//printf("XX superbasic out\n");
2494	}
2495	dualOut_ = dj_[sequenceOut_];
2496	double checkValue = `1.0e-2`;
2497	if (largestDualError_ > `1.0e-5`)
2498	checkValue = `1.0e-1`;
2499	// if stable replace in basis
2500
2501	int updateStatus = factorization_->replaceColumn(this,
2502	rowArray_[`2`],
2503	rowArray_[`1`],
2504	pivotRow_,
2505	alpha_);
2506
2507	// if no pivots, bad update but reasonable alpha - take and invert
2508	if (updateStatus == `2` &&
2509	lastGoodIteration_ == numberIterations_ && fabs(alpha_) > `1.0e-5`)
2510	updateStatus = `4`;
2511	if (updateStatus == `1` \|\| updateStatus == `4`) {
2512	// slight error
2513	if (factorization_->pivots() > `5` \|\| updateStatus == `4`) {
2514	returnCode = -`3`;
2515	}
2516	} else if (updateStatus == `2`) {
2517	// major error
2518	// better to have small tolerance even if slower
2519	factorization_->zeroTolerance(CoinMin(factorization_->zeroTolerance(), `1.0e-15`));
2520	int maxFactor = factorization_->maximumPivots();
2521	if (maxFactor > `10`) {
2522	if (forceFactorization_ < `0`)
2523	forceFactorization_ = maxFactor;
2524	forceFactorization_ = CoinMax(`1`, (forceFactorization_ >> `1`));
2525	}
2526	// later we may need to unwind more e.g. fake bounds
2527	if(lastGoodIteration_ != numberIterations_) {
2528	clearAll();
2529	pivotRow_ = -`1`;
2530	returnCode = -`4`;
2531	} else {
2532	// need to reject something
2533	char x = isColumn(sequenceIn_) ? `'C'` : `'R'`;
2534	handler_->message(CLP_SIMPLEX_FLAG, messages_)
2535	<< x << sequenceWithin(sequenceIn_)
2536	<< CoinMessageEol;
2537	setFlagged(sequenceIn_);
2538	progress_.clearBadTimes();
2539	lastBadIteration_ = numberIterations_; // say be more cautious
2540	clearAll();
2541	pivotRow_ = -`1`;
2542	sequenceOut_ = -`1`;
2543	returnCode = -`5`;
2544
2545	}
2546	return returnCode;
2547	} else if (updateStatus == `3`) {
2548	// out of memory
2549	// increase space if not many iterations
2550	if (factorization_->pivots() <
2551	`0.5` * factorization_->maximumPivots() &&
2552	factorization_->pivots() < `200`)
2553	factorization_->areaFactor(
2554	factorization_->areaFactor() * `1.1`);
2555	returnCode = -`2`; // factorize now
2556	} else if (updateStatus == `5`) {
2557	problemStatus_ = -`2`; // factorize now
2558	}
2559
2560	// update primal solution
2561
2562	double objectiveChange = `0.0`;
2563	// after this rowArray_[1] is not empty - used to update djs
2564	// If pivot row >= numberRows then may be gub
2565	updatePrimalsInPrimal(rowArray_[`1`], theta_, objectiveChange, `1`);
2566
2567	double oldValue = valueIn_;
2568	if (directionIn_ == -`1`) {
2569	// as if from upper bound
2570	if (sequenceIn_ != sequenceOut_) {
2571	// variable becoming basic
2572	valueIn_ -= fabs(theta_);
2573	} else {
2574	valueIn_ = lowerIn_;
2575	}
2576	} else {
2577	// as if from lower bound
2578	if (sequenceIn_ != sequenceOut_) {
2579	// variable becoming basic
2580	valueIn_ += fabs(theta_);
2581	} else {
2582	valueIn_ = upperIn_;
2583	}
2584	}
2585	objectiveChange += dualIn_ * (valueIn_ - oldValue);
2586	// outgoing
2587	if (sequenceIn_ != sequenceOut_) {
2588	if (directionOut_ > `0`) {
2589	valueOut_ = lowerOut_;
2590	} else {
2591	valueOut_ = upperOut_;
2592	}
2593	if(valueOut_ < lower_[sequenceOut_] - primalTolerance_)
2594	valueOut_ = lower_[sequenceOut_] - `0.9` * primalTolerance_;
2595	else if (valueOut_ > upper_[sequenceOut_] + primalTolerance_)
2596	valueOut_ = upper_[sequenceOut_] + `0.9` * primalTolerance_;
2597	// may not be exactly at bound and bounds may have changed
2598	// Make sure outgoing looks feasible
2599	if (!isSuperBasic)
2600	directionOut_ = nonLinearCost_->setOneOutgoing(sequenceOut_, valueOut_);
2601	solution_[sequenceOut_] = valueOut_;
2602	}
2603	// change cost and bounds on incoming if primal
2604	nonLinearCost_->setOne(sequenceIn_, valueIn_);
2605	int whatNext = housekeeping(objectiveChange);
2606	if (keepValue)
2607	solution_[sequenceOut_] = saveValue;
2608	if (isSuperBasic)
2609	setStatus(sequenceOut_, superBasic);
2610	//#define CLP_DEBUG
2611	#if CLP_DEBUG > 1
2612	{
2613	int ninf = matrix_->checkFeasible(this);
2614	if (ninf)
2615	printf("infeas %d\n", ninf);
2616	}
2617	#endif
2618	if (whatNext == `1`) {
2619	returnCode = -`2`; // refactorize
2620	} else if (whatNext == `2`) {
2621	// maximum iterations or equivalent
2622	returnCode = `3`;
2623	} else if(numberIterations_ == lastGoodIteration_
2624	+ `2` * factorization_->maximumPivots()) {
2625	// done a lot of flips - be safe
2626	returnCode = -`2`; // refactorize
2627	}
2628	// Check event
2629	{
2630	int status = eventHandler_->event(ClpEventHandler::endOfIteration);
2631	if (status >= `0`) {
2632	problemStatus_ = `5`;
2633	secondaryStatus_ = ClpEventHandler::endOfIteration;
2634	returnCode = `4`;
2635	}
2636	}
2637	return returnCode;
2638	}
2639	// A sequential LP method
2640	int
2641	ClpSimplexNonlinear::primalSLP(int numberPasses, double deltaTolerance)
2642	{
2643	// Are we minimizing or maximizing
2644	double whichWay = optimizationDirection();
2645	if (whichWay < `0.0`)
2646	whichWay = -`1.0`;
2647	else if (whichWay > `0.0`)
2648	whichWay = `1.0`;
2649
2650
2651	int numberColumns = this->numberColumns();
2652	int numberRows = this->numberRows();
2653	double * columnLower = this->columnLower();
2654	double * columnUpper = this->columnUpper();
2655	double * solution = this->primalColumnSolution();
2656
2657	if (objective_->type() < `2`) {
2658	// no nonlinear part
2659	return ClpSimplex::primal(`0`);
2660	}
2661	// Get list of non linear columns
2662	char * markNonlinear = new char[numberColumns];
2663	memset(markNonlinear, `0`, numberColumns);
2664	int numberNonLinearColumns = objective_->markNonlinear(markNonlinear);
2665
2666	if (!numberNonLinearColumns) {
2667	delete [] markNonlinear;
2668	// no nonlinear part
2669	return ClpSimplex::primal(`0`);
2670	}
2671	int iColumn;
2672	int * listNonLinearColumn = new int[numberNonLinearColumns];
2673	numberNonLinearColumns = `0`;
2674	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
2675	if(markNonlinear[iColumn])
2676	listNonLinearColumn[numberNonLinearColumns++] = iColumn;
2677	}
2678	// Replace objective
2679	ClpObjective * trueObjective = objective_;
2680	objective_ = new ClpLinearObjective (NULL, numberColumns);
2681	double * objective = this->objective();
2682
2683	// get feasible
2684	if (this->status() < `0` \|\| numberPrimalInfeasibilities())
2685	ClpSimplex::primal(`1`);
2686	// still infeasible
2687	if (numberPrimalInfeasibilities()) {
2688	delete [] listNonLinearColumn;
2689	return `0`;
2690	}
2691	algorithm_ = `1`;
2692	int jNon;
2693	int * last[`3`];
2694
2695	double * trust = new double[numberNonLinearColumns];
2696	double * trueLower = new double[numberNonLinearColumns];
2697	double * trueUpper = new double[numberNonLinearColumns];
2698	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
2699	iColumn = listNonLinearColumn[jNon];
2700	trust[jNon] = `0.5`;
2701	trueLower[jNon] = columnLower[iColumn];
2702	trueUpper[jNon] = columnUpper[iColumn];
2703	if (solution[iColumn] < trueLower[jNon])
2704	solution[iColumn] = trueLower[jNon];
2705	else if (solution[iColumn] > trueUpper[jNon])
2706	solution[iColumn] = trueUpper[jNon];
2707	}
2708	int saveLogLevel = logLevel();
2709	int iPass;
2710	double lastObjective = `1.0e31`;
2711	double * saveSolution = new double [numberColumns];
2712	double * saveRowSolution = new double [numberRows];
2713	memset(saveRowSolution, `0`, numberRows * sizeof(double));
2714	double * savePi = new double [numberRows];
2715	double * safeSolution = new double [numberColumns];
2716	unsigned char * saveStatus = new unsigned char[numberRows+numberColumns];
2717	#define MULTIPLE 0
2718	#if MULTIPLE>2
2719	// Duplication but doesn't really matter
2720	double * saveSolutionM[MULTIPLE
2721	};
2722	for (jNon=`0`; jNon<MULTIPLE; jNon++)
2723	{
2724	saveSolutionM[jNon] = new double[numberColumns];
2725	CoinMemcpyN(solution, numberColumns, saveSolutionM);
2726	}
2727	#endif
2728	double targetDrop = `1.0e31`;
2729	double objectiveOffset;
2730	getDblParam(ClpObjOffset, objectiveOffset);
2731	// 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
2732	for (iPass = `0`; iPass < `3`; iPass++)
2733	{
2734	last[iPass] = new int[numberNonLinearColumns];
2735	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++)
2736	last[iPass][jNon] = `0`;
2737	}
2738	// goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
2739	int goodMove = -`2`;
2740	char * statusCheck = new char[numberColumns];
2741	double * changeRegion = new double [numberColumns];
2742	double offset = `0.0`;
2743	double objValue = `0.0`;
2744	int exitPass = `2` * numberPasses + `10`;
2745	for (iPass = `0`; iPass < numberPasses; iPass++)
2746	{
2747	exitPass--;
2748	// redo objective
2749	offset = `0.0`;
2750	objValue = -objectiveOffset;
2751	// make sure x updated
2752	trueObjective->newXValues();
2753	double theta = -`1.0`;
2754	double maxTheta = COIN_DBL_MAX;
2755	//maxTheta=1.0;
2756	if (iPass) {
2757	int jNon = `0`;
2758	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
2759	changeRegion[iColumn] = solution[iColumn] - saveSolution[iColumn];
2760	double alpha = changeRegion[iColumn];
2761	double oldValue = saveSolution[iColumn];
2762	if (markNonlinear[iColumn] == `0`) {
2763	// linear
2764	if (alpha < -`1.0e-15`) {
2765	// variable going towards lower bound
2766	double bound = columnLower[iColumn];
2767	oldValue -= bound;
2768	if (oldValue + maxTheta * alpha < `0.0`) {
2769	maxTheta = CoinMax(`0.0`, oldValue / (-alpha));
2770	}
2771	} else if (alpha > `1.0e-15`) {
2772	// variable going towards upper bound
2773	double bound = columnUpper[iColumn];
2774	oldValue = bound - oldValue;
2775	if (oldValue - maxTheta * alpha < `0.0`) {
2776	maxTheta = CoinMax(`0.0`, oldValue / alpha);
2777	}
2778	}
2779	} else {
2780	// nonlinear
2781	if (alpha < -`1.0e-15`) {
2782	// variable going towards lower bound
2783	double bound = trueLower[jNon];
2784	oldValue -= bound;
2785	if (oldValue + maxTheta * alpha < `0.0`) {
2786	maxTheta = CoinMax(`0.0`, oldValue / (-alpha));
2787	}
2788	} else if (alpha > `1.0e-15`) {
2789	// variable going towards upper bound
2790	double bound = trueUpper[jNon];
2791	oldValue = bound - oldValue;
2792	if (oldValue - maxTheta * alpha < `0.0`) {
2793	maxTheta = CoinMax(`0.0`, oldValue / alpha);
2794	}
2795	}
2796	jNon++;
2797	}
2798	}
2799	// make sure both accurate
2800	memset(rowActivity_, `0`, numberRows_ * sizeof(double));
2801	times(`1.0`, solution, rowActivity_);
2802	memset(saveRowSolution, `0`, numberRows_ * sizeof(double));
2803	times(`1.0`, saveSolution, saveRowSolution);
2804	for (int iRow = `0`; iRow < numberRows; iRow++) {
2805	double alpha = rowActivity_[iRow] - saveRowSolution[iRow];
2806	double oldValue = saveRowSolution[iRow];
2807	if (alpha < -`1.0e-15`) {
2808	// variable going towards lower bound
2809	double bound = rowLower_[iRow];
2810	oldValue -= bound;
2811	if (oldValue + maxTheta * alpha < `0.0`) {
2812	maxTheta = CoinMax(`0.0`, oldValue / (-alpha));
2813	}
2814	} else if (alpha > `1.0e-15`) {
2815	// variable going towards upper bound
2816	double bound = rowUpper_[iRow];
2817	oldValue = bound - oldValue;
2818	if (oldValue - maxTheta * alpha < `0.0`) {
2819	maxTheta = CoinMax(`0.0`, oldValue / alpha);
2820	}
2821	}
2822	}
2823	} else {
2824	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
2825	changeRegion[iColumn] = `0.0`;
2826	saveSolution[iColumn] = solution[iColumn];
2827	}
2828	CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
2829	}
2830	// get current value anyway
2831	double predictedObj, thetaObj;
2832	double maxTheta2 = `2.0`; // to work out a b c
2833	double theta2 = trueObjective->stepLength(this, saveSolution, changeRegion, maxTheta2,
2834	objValue, predictedObj, thetaObj);
2835	int lastMoveStatus = goodMove;
2836	if (goodMove >= `0`) {
2837	theta = CoinMin(theta2, maxTheta);
2838	#ifdef CLP_DEBUG
2839	if (handler_->logLevel() & `32`)
2840	printf("theta %g, current %g, at maxtheta %g, predicted %g\n",
2841	theta, objValue, thetaObj, predictedObj);
2842	#endif
2843	if (theta > `0.0` && theta <= `1.0`) {
2844	// update solution
2845	double lambda = `1.0` - theta;
2846	for (iColumn = `0`; iColumn < numberColumns; iColumn++)
2847	solution[iColumn] = lambda * saveSolution[iColumn]
2848	+ theta * solution[iColumn];
2849	memset(rowActivity_, `0`, numberRows_ * sizeof(double));
2850	times(`1.0`, solution, rowActivity_);
2851	if (lambda > `0.999`) {
2852	CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
2853	CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
2854	}
2855	// Do local minimization
2856	#define LOCAL
2857	#ifdef LOCAL
2858	bool absolutelyOptimal = false;
2859	int saveScaling = scalingFlag();
2860	scaling(`0`);
2861	int savePerturbation = perturbation_;
2862	perturbation_ = `100`;
2863	if (saveLogLevel == `1`)
2864	setLogLevel(`0`);
2865	int status = startup(`1`);
2866	if (!status) {
2867	int numberTotal = numberRows_ + numberColumns_;
2868	// resize arrays
2869	for (int i = `0`; i < `4`; i++) {
2870	rowArray_[i]->reserve(CoinMax(numberRows_ + numberColumns_, rowArray_[i]->capacity()));
2871	}
2872	CoinIndexedVector * longArray = rowArray_[`3`];
2873	CoinIndexedVector * rowArray = rowArray_[`0`];
2874	//CoinIndexedVector columnArray = columnArray_[0];*
2875	CoinIndexedVector * spare = rowArray_[`1`];
2876	double * work = longArray->denseVector();
2877	int * which = longArray->getIndices();
2878	int nPass = `100`;
2879	//bool conjugate=false;
2880	// Put back true bounds
2881	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
2882	int iColumn = listNonLinearColumn[jNon];
2883	double value;
2884	value = trueLower[jNon];
2885	trueLower[jNon] = lower_[iColumn];
2886	lower_[iColumn] = value;
2887	value = trueUpper[jNon];
2888	trueUpper[jNon] = upper_[iColumn];
2889	upper_[iColumn] = value;
2890	switch(getStatus(iColumn)) {
2891
2892	case basic:
2893	case isFree:
2894	case superBasic:
2895	break;
2896	case isFixed:
2897	case atUpperBound:
2898	case atLowerBound:
2899	if (solution_[iColumn] > lower_[iColumn] + primalTolerance_) {
2900	if(solution_[iColumn] < upper_[iColumn] - primalTolerance_) {
2901	setStatus(iColumn, superBasic);
2902	} else {
2903	setStatus(iColumn, atUpperBound);
2904	}
2905	} else {
2906	if(solution_[iColumn] < upper_[iColumn] - primalTolerance_) {
2907	setStatus(iColumn, atLowerBound);
2908	} else {
2909	setStatus(iColumn, isFixed);
2910	}
2911	}
2912	break;
2913	}
2914	}
2915	for (int jPass = `0`; jPass < nPass; jPass++) {
2916	trueObjective->reducedGradient(this, dj_, true);
2917	// get direction vector in longArray
2918	longArray->clear();
2919	double normFlagged = `0.0`;
2920	double normUnflagged = `0.0`;
2921	int numberNonBasic = `0`;
2922	directionVector(longArray, spare, rowArray, `0`,
2923	normFlagged, normUnflagged, numberNonBasic);
2924	if (normFlagged + normUnflagged < `1.0e-8`) {
2925	absolutelyOptimal = true;
2926	break; //looks optimal
2927	}
2928	double djNorm = `0.0`;
2929	int iIndex;
2930	for (iIndex = `0`; iIndex < numberNonBasic; iIndex++) {
2931	int iSequence = which[iIndex];
2932	double alpha = work[iSequence];
2933	djNorm += alpha * alpha;
2934	}
2935	//if (!jPass)
2936	//printf("dj norm %g - %d \n",djNorm,numberNonBasic);
2937	//int number=longArray->getNumElements();
2938	if (!numberNonBasic) {
2939	// looks optimal
2940	absolutelyOptimal = true;
2941	break;
2942	}
2943	int bestSequence = -`1`;
2944	double theta = `1.0e30`;
2945	int iSequence;
2946	for (iSequence = `0`; iSequence < numberTotal; iSequence++) {
2947	double alpha = work[iSequence];
2948	double oldValue = solution_[iSequence];
2949	if (alpha < -`1.0e-15`) {
2950	// variable going towards lower bound
2951	double bound = lower_[iSequence];
2952	oldValue -= bound;
2953	if (oldValue + theta * alpha < `0.0`) {
2954	bestSequence = iSequence;
2955	theta = CoinMax(`0.0`, oldValue / (-alpha));
2956	}
2957	} else if (alpha > `1.0e-15`) {
2958	// variable going towards upper bound
2959	double bound = upper_[iSequence];
2960	oldValue = bound - oldValue;
2961	if (oldValue - theta * alpha < `0.0`) {
2962	bestSequence = iSequence;
2963	theta = CoinMax(`0.0`, oldValue / alpha);
2964	}
2965	}
2966	}
2967	// Now find minimum of function
2968	double currentObj;
2969	double predictedObj;
2970	double thetaObj;
2971	// need to use true objective
2972	double * saveCost = cost_;
2973	cost_ = NULL;
2974	double objTheta = trueObjective->stepLength(this, solution_, work, theta,
2975	currentObj, predictedObj, thetaObj);
2976	cost_ = saveCost;
2977	#ifdef CLP_DEBUG
2978	if (handler_->logLevel() & `32`)
2979	printf("current obj %g thetaObj %g, predictedObj %g\n", currentObj, thetaObj, predictedObj);
2980	#endif
2981	//printf("current obj %g thetaObj %g, predictedObj %g\n",currentObj,thetaObj,predictedObj);
2982	//printf("objTheta %g theta %g\n",objTheta,theta);
2983	if (theta > objTheta) {
2984	theta = objTheta;
2985	thetaObj = predictedObj;
2986	}
2987	// update one used outside
2988	objValue = currentObj;
2989	if (theta > `1.0e-9` &&
2990	(currentObj - thetaObj < -CoinMax(`1.0e-8`, `1.0e-15` * fabs(currentObj)) \|\| jPass < `5`)) {
2991	// Update solution
2992	for (iSequence = `0`; iSequence < numberTotal; iSequence++) {
2993	double alpha = work[iSequence];
2994	if (alpha) {
2995	double value = solution_[iSequence] + theta * alpha;
2996	solution_[iSequence] = value;
2997	switch(getStatus(iSequence)) {
2998
2999	case basic:
3000	case isFixed:
3001	case isFree:
3002	break;
3003	case atUpperBound:
3004	case atLowerBound:
3005	case superBasic:
3006	if (value > lower_[iSequence] + primalTolerance_) {
3007	if(value < upper_[iSequence] - primalTolerance_) {
3008	setStatus(iSequence, superBasic);
3009	} else {
3010	setStatus(iSequence, atUpperBound);
3011	}
3012	} else {
3013	if(value < upper_[iSequence] - primalTolerance_) {
3014	setStatus(iSequence, atLowerBound);
3015	} else {
3016	setStatus(iSequence, isFixed);
3017	}
3018	}
3019	break;
3020	}
3021	}
3022	}
3023	} else {
3024	break;
3025	}
3026	}
3027	// Put back fake bounds
3028	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3029	int iColumn = listNonLinearColumn[jNon];
3030	double value;
3031	value = trueLower[jNon];
3032	trueLower[jNon] = lower_[iColumn];
3033	lower_[iColumn] = value;
3034	value = trueUpper[jNon];
3035	trueUpper[jNon] = upper_[iColumn];
3036	upper_[iColumn] = value;
3037	}
3038	}
3039	problemStatus_ = `0`;
3040	finish();
3041	scaling(saveScaling);
3042	perturbation_ = savePerturbation;
3043	setLogLevel(saveLogLevel);
3044	#endif
3045	// redo rowActivity
3046	memset(rowActivity_, `0`, numberRows_ * sizeof(double));
3047	times(`1.0`, solution, rowActivity_);
3048	if (theta > `0.99999` && theta2 < `1.9` && !absolutelyOptimal) {
3049	// If big changes then tighten
3050	/ thetaObj is objvalue + a22 +b2
3051	predictedObj is objvalue + atheta2theta2 +btheta2*
3052	*/
3053	double rhs1 = thetaObj - objValue;
3054	double rhs2 = predictedObj - objValue;
3055	double subtractB = theta2 * `0.5`;
3056	double a = (rhs2 - subtractB * rhs1) / (theta2 * theta2 - `4.0` * subtractB);
3057	double b = `0.5` * (rhs1 - `4.0` * a);
3058	if (fabs(a + b) > `1.0e-2`) {
3059	// tighten all
3060	goodMove = -`1`;
3061	}
3062	}
3063	}
3064	}
3065	CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, `2`), numberColumns,
3066	objective);
3067	//printf("offset comp %g orig %g\n",offset,objectiveOffset);
3068	// could do faster
3069	trueObjective->stepLength(this, solution, changeRegion, `0.0`,
3070	objValue, predictedObj, thetaObj);
3071	#ifdef CLP_INVESTIGATE
3072	printf("offset comp %g orig %g - obj from stepLength %g\n", offset, objectiveOffset, objValue);
3073	#endif
3074	setDblParam(ClpObjOffset, objectiveOffset + offset);
3075	int * temp = last[`2`];
3076	last[`2`] = last[`1`];
3077	last[`1`] = last[`0`];
3078	last[`0`] = temp;
3079	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3080	iColumn = listNonLinearColumn[jNon];
3081	double change = solution[iColumn] - saveSolution[iColumn];
3082	if (change < -`1.0e-5`) {
3083	if (fabs(change + trust[jNon]) < `1.0e-5`)
3084	temp[jNon] = -`1`;
3085	else
3086	temp[jNon] = -`2`;
3087	} else if(change > `1.0e-5`) {
3088	if (fabs(change - trust[jNon]) < `1.0e-5`)
3089	temp[jNon] = `1`;
3090	else
3091	temp[jNon] = `2`;
3092	} else {
3093	temp[jNon] = `0`;
3094	}
3095	}
3096	// goodMove +1 yes, 0 no, -1 last was bad - just halve gaps, -2 do nothing
3097	double maxDelta = `0.0`;
3098	if (goodMove >= `0`) {
3099	if (objValue - lastObjective <= `1.0e-15` * fabs(lastObjective))
3100	goodMove = `1`;
3101	else
3102	goodMove = `0`;
3103	} else {
3104	maxDelta = `1.0e10`;
3105	}
3106	double maxGap = `0.0`;
3107	int numberSmaller = `0`;
3108	int numberSmaller2 = `0`;
3109	int numberLarger = `0`;
3110	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3111	iColumn = listNonLinearColumn[jNon];
3112	maxDelta = CoinMax(maxDelta,
3113	fabs(solution[iColumn] - saveSolution[iColumn]));
3114	if (goodMove > `0`) {
3115	if (last[`0`][jNon]*last[`1`][jNon] < `0`) {
3116	// halve
3117	trust[jNon] *= `0.5`;
3118	numberSmaller2++;
3119	} else {
3120	if (last[`0`][jNon] == last[`1`][jNon] &&
3121	last[`0`][jNon] == last[`2`][jNon])
3122	trust[jNon] = CoinMin(`1.5` * trust[jNon], `1.0e6`);
3123	numberLarger++;
3124	}
3125	} else if (goodMove != -`2` && trust[jNon] > `10.0` * deltaTolerance) {
3126	trust[jNon] *= `0.2`;
3127	numberSmaller++;
3128	}
3129	maxGap = CoinMax(maxGap, trust[jNon]);
3130	}
3131	#ifdef CLP_DEBUG
3132	if (handler_->logLevel() & `32`)
3133	std::cout << "largest gap is " << maxGap << " "
3134	<< numberSmaller + numberSmaller2 << " reduced ("
3135	<< numberSmaller << " badMove ), "
3136	<< numberLarger << " increased" << std::endl;
3137	#endif
3138	if (iPass > `10000`) {
3139	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++)
3140	trust[jNon] *= `0.0001`;
3141	}
3142	if(lastMoveStatus == -`1` && goodMove == -`1`)
3143	goodMove = `1`; // to force solve
3144	if (goodMove > `0`) {
3145	double drop = lastObjective - objValue;
3146	handler_->message(CLP_SLP_ITER, messages_)
3147	<< iPass << objValue - objectiveOffset
3148	<< drop << maxDelta
3149	<< CoinMessageEol;
3150	if (iPass > `20` && drop < `1.0e-12` * fabs(objValue))
3151	drop = `0.999e-4`; // so will exit
3152	if (maxDelta < deltaTolerance && drop < `1.0e-4` && goodMove && theta < `0.99999`) {
3153	if (handler_->logLevel() > `1`)
3154	std::cout << "Exiting as maxDelta < tolerance and small drop" << std::endl;
3155	break;
3156	}
3157	} else if (!numberSmaller && iPass > `1`) {
3158	if (handler_->logLevel() > `1`)
3159	std::cout << "Exiting as all gaps small" << std::endl;
3160	break;
3161	}
3162	if (!iPass)
3163	goodMove = `1`;
3164	targetDrop = `0.0`;
3165	double * r = this->dualColumnSolution();
3166	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3167	iColumn = listNonLinearColumn[jNon];
3168	columnLower[iColumn] = CoinMax(solution[iColumn]
3169	- trust[jNon],
3170	trueLower[jNon]);
3171	columnUpper[iColumn] = CoinMin(solution[iColumn]
3172	+ trust[jNon],
3173	trueUpper[jNon]);
3174	}
3175	if (iPass) {
3176	// get reduced costs
3177	this->matrix()->transposeTimes(savePi,
3178	this->dualColumnSolution());
3179	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3180	iColumn = listNonLinearColumn[jNon];
3181	double dj = objective[iColumn] - r[iColumn];
3182	r[iColumn] = dj;
3183	if (dj < -dualTolerance_)
3184	targetDrop -= dj * (columnUpper[iColumn] - solution[iColumn]);
3185	else if (dj > dualTolerance_)
3186	targetDrop -= dj * (columnLower[iColumn] - solution[iColumn]);
3187	}
3188	} else {
3189	memset(r, `0`, numberColumns * sizeof(double));
3190	}
3191	#if 0
3192	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3193	iColumn = listNonLinearColumn[jNon];
3194	if (statusCheck[iColumn] == `'L'` && r[iColumn] < -`1.0e-4`) {
3195	columnLower[iColumn] = CoinMax(solution[iColumn],
3196	trueLower[jNon]);
3197	columnUpper[iColumn] = CoinMin(solution[iColumn]
3198	+ trust[jNon],
3199	trueUpper[jNon]);
3200	} else if (statusCheck[iColumn] == `'U'` && r[iColumn] > `1.0e-4`) {
3201	columnLower[iColumn] = CoinMax(solution[iColumn]
3202	- trust[jNon],
3203	trueLower[jNon]);
3204	columnUpper[iColumn] = CoinMin(solution[iColumn],
3205	trueUpper[jNon]);
3206	} else {
3207	columnLower[iColumn] = CoinMax(solution[iColumn]
3208	- trust[jNon],
3209	trueLower[jNon]);
3210	columnUpper[iColumn] = CoinMin(solution[iColumn]
3211	+ trust[jNon],
3212	trueUpper[jNon]);
3213	}
3214	}
3215	#endif
3216	if (goodMove > `0`) {
3217	CoinMemcpyN(solution, numberColumns, saveSolution);
3218	CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
3219	CoinMemcpyN(this->dualRowSolution(), numberRows, savePi);
3220	CoinMemcpyN(status_, numberRows + numberColumns, saveStatus);
3221	#if MULTIPLE>2
3222	double * tempSol = saveSolutionM[`0`];
3223	for (jNon = `0`; jNon < MULTIPLE - `1`; jNon++) {
3224	saveSolutionM[jNon] = saveSolutionM[jNon+`1`];
3225	}
3226	saveSolutionM[MULTIPLE-`1`] = tempSol;
3227	CoinMemcpyN(solution, numberColumns, tempSol);
3228
3229	#endif
3230
3231	#ifdef CLP_DEBUG
3232	if (handler_->logLevel() & `32`)
3233	std::cout << "Pass - " << iPass
3234	<< ", target drop is " << targetDrop
3235	<< std::endl;
3236	#endif
3237	lastObjective = objValue;
3238	if (targetDrop < CoinMax(`1.0e-8`, CoinMin(`1.0e-6`, `1.0e-6` * fabs(objValue))) && goodMove && iPass > `3`) {
3239	if (handler_->logLevel() > `1`)
3240	printf("Exiting on target drop %g\n", targetDrop);
3241	break;
3242	}
3243	#ifdef CLP_DEBUG
3244	{
3245	double * r = this->dualColumnSolution();
3246	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3247	iColumn = listNonLinearColumn[jNon];
3248	if (handler_->logLevel() & `32`)
3249	printf("Trust %d %g - solution %d %g obj %g dj %g state %c - bounds %g %g\n",
3250	jNon, trust[jNon], iColumn, solution[iColumn], objective[iColumn],
3251	r[iColumn], statusCheck[iColumn], columnLower[iColumn],
3252	columnUpper[iColumn]);
3253	}
3254	}
3255	#endif
3256	//setLogLevel(63);
3257	this->scaling(false);
3258	if (saveLogLevel == `1`)
3259	setLogLevel(`0`);
3260	ClpSimplex::primal(`1`);
3261	algorithm_ = `1`;
3262	setLogLevel(saveLogLevel);
3263	#ifdef CLP_DEBUG
3264	if (this->status()) {
3265	writeMps("xx.mps");
3266	}
3267	#endif
3268	if (this->status() == `1`) {
3269	// not feasible ! - backtrack and exit
3270	// use safe solution
3271	CoinMemcpyN(safeSolution, numberColumns, solution);
3272	CoinMemcpyN(solution, numberColumns, saveSolution);
3273	memset(rowActivity_, `0`, numberRows_ * sizeof(double));
3274	times(`1.0`, solution, rowActivity_);
3275	CoinMemcpyN(rowActivity_, numberRows, saveRowSolution);
3276	CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
3277	CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
3278	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3279	iColumn = listNonLinearColumn[jNon];
3280	columnLower[iColumn] = CoinMax(solution[iColumn]
3281	- trust[jNon],
3282	trueLower[jNon]);
3283	columnUpper[iColumn] = CoinMin(solution[iColumn]
3284	+ trust[jNon],
3285	trueUpper[jNon]);
3286	}
3287	break;
3288	} else {
3289	// save in case problems
3290	CoinMemcpyN(solution, numberColumns, safeSolution);
3291	}
3292	if (problemStatus_ == `3`)
3293	break; // probably user interrupt
3294	goodMove = `1`;
3295	} else {
3296	// bad pass - restore solution
3297	#ifdef CLP_DEBUG
3298	if (handler_->logLevel() & `32`)
3299	printf("Backtracking\n");
3300	#endif
3301	CoinMemcpyN(saveSolution, numberColumns, solution);
3302	CoinMemcpyN(saveRowSolution, numberRows, rowActivity_);
3303	CoinMemcpyN(savePi, numberRows, this->dualRowSolution());
3304	CoinMemcpyN(saveStatus, numberRows + numberColumns, status_);
3305	iPass--;
3306	assert (exitPass > `0`);
3307	goodMove = -`1`;
3308	}
3309	}
3310	#if MULTIPLE>2
3311	for (jNon = `0`; jNon < MULTIPLE; jNon++)
3312	{
3313	delete [] saveSolutionM[jNon];
3314	}
3315	#endif
3316	// restore solution
3317	CoinMemcpyN(saveSolution, numberColumns, solution);
3318	CoinMemcpyN(saveRowSolution, numberRows, rowActivity_);
3319	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++)
3320	{
3321	iColumn = listNonLinearColumn[jNon];
3322	columnLower[iColumn] = CoinMax(solution[iColumn],
3323	trueLower[jNon]);
3324	columnUpper[iColumn] = CoinMin(solution[iColumn],
3325	trueUpper[jNon]);
3326	}
3327	delete [] markNonlinear;
3328	delete [] statusCheck;
3329	delete [] savePi;
3330	delete [] saveStatus;
3331	// redo objective
3332	CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, `2`), numberColumns,
3333	objective);
3334	ClpSimplex::primal(`1`);
3335	delete objective_;
3336	objective_ = trueObjective;
3337	// redo values
3338	setDblParam(ClpObjOffset, objectiveOffset);
3339	objectiveValue_ += whichWay * offset;
3340	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++)
3341	{
3342	iColumn = listNonLinearColumn[jNon];
3343	columnLower[iColumn] = trueLower[jNon];
3344	columnUpper[iColumn] = trueUpper[jNon];
3345	}
3346	delete [] saveSolution;
3347	delete [] safeSolution;
3348	delete [] saveRowSolution;
3349	for (iPass = `0`; iPass < `3`; iPass++)
3350	delete [] last[iPass];
3351	delete [] trust;
3352	delete [] trueUpper;
3353	delete [] trueLower;
3354	delete [] listNonLinearColumn;
3355	delete [] changeRegion;
3356	// temp
3357	//setLogLevel(63);
3358	return `0`;
3359	}
3360	/ Primal algorithm for nonlinear constraints*
3361	Using a semi-trust region approach as for pooling problem
3362	This is in because I have it lying around
3363
3364	*/
3365	int
3366	ClpSimplexNonlinear::primalSLP(int numberConstraints, ClpConstraint ** constraints,
3367	int numberPasses, double deltaTolerance)
3368	{
3369	if (!numberConstraints) {
3370	// no nonlinear constraints - may be nonlinear objective
3371	return primalSLP(numberPasses, deltaTolerance);
3372	}
3373	// Are we minimizing or maximizing
3374	double whichWay = optimizationDirection();
3375	if (whichWay < `0.0`)
3376	whichWay = -`1.0`;
3377	else if (whichWay > `0.0`)
3378	whichWay = `1.0`;
3379	// check all matrix for odd rows is empty
3380	int iConstraint;
3381	int numberBad = `0`;
3382	CoinPackedMatrix * columnCopy = matrix();
3383	// Get a row copy in standard format
3384	CoinPackedMatrix copy;
3385	copy.reverseOrderedCopyOf(*columnCopy);
3386	// get matrix data pointers
3387	//const int column = copy.getIndices();*
3388	//const CoinBigIndex rowStart = copy.getVectorStarts();*
3389	const int * rowLength = copy.getVectorLengths();
3390	//const double elementByRow = copy.getElements();*
3391	int numberArtificials = `0`;
3392	// We could use nonlinearcost to do segments - maybe later
3393	#define SEGMENTS 3
3394	// Penalties may be adjusted by duals
3395	// Both these should be modified depending on problem
3396	// Possibly start with big bounds
3397	//double penalties[]={1.0e-3,1.0e7,1.0e9};
3398	double penalties[] = {`1.0e7`, `1.0e8`, `1.0e9`};
3399	double bounds[] = {`1.0e-2`, `1.0e2`, COIN_DBL_MAX};
3400	// see how many extra we need
3401	CoinBigIndex numberExtra = `0`;
3402	for (iConstraint = `0`; iConstraint < numberConstraints; iConstraint++) {
3403	ClpConstraint * constraint = constraints[iConstraint];
3404	int iRow = constraint->rowNumber();
3405	assert (iRow >= `0`);
3406	int n = constraint->numberCoefficients() - rowLength[iRow];
3407	numberExtra += n;
3408	if (iRow >= numberRows_)
3409	numberBad++;
3410	else if (rowLength[iRow] && n)
3411	numberBad++;
3412	if (rowLower_[iRow] > -`1.0e20`)
3413	numberArtificials++;
3414	if (rowUpper_[iRow] < `1.0e20`)
3415	numberArtificials++;
3416	}
3417	if (numberBad)
3418	return numberBad;
3419	ClpObjective * trueObjective = NULL;
3420	if (objective_->type() >= `2`) {
3421	// Replace objective
3422	trueObjective = objective_;
3423	objective_ = new ClpLinearObjective (NULL, numberColumns_);
3424	}
3425	ClpSimplex newModel(*this);
3426	// we can put back true objective
3427	if (trueObjective) {
3428	// Replace objective
3429	delete objective_;
3430	objective_ = trueObjective;
3431	}
3432	int numberColumns2 = numberColumns_;
3433	int numberSmallGap = numberArtificials;
3434	if (numberArtificials) {
3435	numberArtificials *= SEGMENTS;
3436	numberColumns2 += numberArtificials;
3437	int * addStarts = new int [numberArtificials+`1`];
3438	int * addRow = new int[numberArtificials];
3439	double * addElement = new double[numberArtificials];
3440	double * addUpper = new double[numberArtificials];
3441	addStarts[`0`] = `0`;
3442	double * addCost = new double [numberArtificials];
3443	numberArtificials = `0`;
3444	for (iConstraint = `0`; iConstraint < numberConstraints; iConstraint++) {
3445	ClpConstraint * constraint = constraints[iConstraint];
3446	int iRow = constraint->rowNumber();
3447	if (rowLower_[iRow] > -`1.0e20`) {
3448	for (int k = `0`; k < SEGMENTS; k++) {
3449	addRow[numberArtificials] = iRow;
3450	addElement[numberArtificials] = `1.0`;
3451	addCost[numberArtificials] = penalties[k];
3452	addUpper[numberArtificials] = bounds[k];
3453	numberArtificials++;
3454	addStarts[numberArtificials] = numberArtificials;
3455	}
3456	}
3457	if (rowUpper_[iRow] < `1.0e20`) {
3458	for (int k = `0`; k < SEGMENTS; k++) {
3459	addRow[numberArtificials] = iRow;
3460	addElement[numberArtificials] = -`1.0`;
3461	addCost[numberArtificials] = penalties[k];
3462	addUpper[numberArtificials] = bounds[k];
3463	numberArtificials++;
3464	addStarts[numberArtificials] = numberArtificials;
3465	}
3466	}
3467	}
3468	newModel.addColumns(numberArtificials, NULL, addUpper, addCost,
3469	addStarts, addRow, addElement);
3470	delete [] addStarts;
3471	delete [] addRow;
3472	delete [] addElement;
3473	delete [] addUpper;
3474	delete [] addCost;
3475	// newModel.primal(1);
3476	}
3477	// find nonlinear columns
3478	int * listNonLinearColumn = new int [numberColumns_+numberSmallGap];
3479	char * mark = new char[numberColumns_];
3480	memset(mark, `0`, numberColumns_);
3481	for (iConstraint = `0`; iConstraint < numberConstraints; iConstraint++) {
3482	ClpConstraint * constraint = constraints[iConstraint];
3483	constraint->markNonlinear(mark);
3484	}
3485	if (trueObjective)
3486	trueObjective->markNonlinear(mark);
3487	int iColumn;
3488	int numberNonLinearColumns = `0`;
3489	for (iColumn = `0`; iColumn < numberColumns_; iColumn++) {
3490	if (mark[iColumn])
3491	listNonLinearColumn[numberNonLinearColumns++] = iColumn;
3492	}
3493	// and small gap artificials
3494	for (iColumn = numberColumns_; iColumn < numberColumns2; iColumn += SEGMENTS) {
3495	listNonLinearColumn[numberNonLinearColumns++] = iColumn;
3496	}
3497	// build row copy of matrix with space for nonzeros
3498	// Get column copy
3499	columnCopy = newModel.matrix();
3500	copy.reverseOrderedCopyOf(*columnCopy);
3501	// get matrix data pointers
3502	const int * column = copy.getIndices();
3503	const CoinBigIndex * rowStart = copy.getVectorStarts();
3504	rowLength = copy.getVectorLengths();
3505	const double * elementByRow = copy.getElements();
3506	numberExtra += copy.getNumElements();
3507	CoinBigIndex * newStarts = new CoinBigIndex [numberRows_+`1`];
3508	int * newColumn = new int[numberExtra];
3509	double * newElement = new double[numberExtra];
3510	newStarts[`0`] = `0`;
3511	int * backRow = new int [numberRows_];
3512	int iRow;
3513	for (iRow = `0`; iRow < numberRows_; iRow++)
3514	backRow[iRow] = -`1`;
3515	for (iConstraint = `0`; iConstraint < numberConstraints; iConstraint++) {
3516	ClpConstraint * constraint = constraints[iConstraint];
3517	iRow = constraint->rowNumber();
3518	backRow[iRow] = iConstraint;
3519	}
3520	numberExtra = `0`;
3521	for (iRow = `0`; iRow < numberRows_; iRow++) {
3522	if (backRow[iRow] < `0`) {
3523	// copy normal
3524	for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow];
3525	j++) {
3526	newColumn[numberExtra] = column[j];
3527	newElement[numberExtra++] = elementByRow[j];
3528	}
3529	} else {
3530	ClpConstraint * constraint = constraints[backRow[iRow]];
3531	assert(iRow == constraint->rowNumber());
3532	int numberArtificials = `0`;
3533	if (rowLower_[iRow] > -`1.0e20`)
3534	numberArtificials += SEGMENTS;
3535	if (rowUpper_[iRow] < `1.0e20`)
3536	numberArtificials += SEGMENTS;
3537	if (numberArtificials == rowLength[iRow]) {
3538	// all possible
3539	memset(mark, `0`, numberColumns_);
3540	constraint->markNonzero(mark);
3541	for (int k = `0`; k < numberColumns_; k++) {
3542	if (mark[k]) {
3543	newColumn[numberExtra] = k;
3544	newElement[numberExtra++] = `1.0`;
3545	}
3546	}
3547	// copy artificials
3548	for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow];
3549	j++) {
3550	newColumn[numberExtra] = column[j];
3551	newElement[numberExtra++] = elementByRow[j];
3552	}
3553	} else {
3554	// there already
3555	// copy
3556	for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow];
3557	j++) {
3558	newColumn[numberExtra] = column[j];
3559	assert (elementByRow[j]);
3560	newElement[numberExtra++] = elementByRow[j];
3561	}
3562	}
3563	}
3564	newStarts[iRow+`1`] = numberExtra;
3565	}
3566	delete [] backRow;
3567	CoinPackedMatrix saveMatrix(false, numberColumns2, numberRows_,
3568	numberExtra, newElement, newColumn, newStarts, NULL, `0.0`, `0.0`);
3569	delete [] newStarts;
3570	delete [] newColumn;
3571	delete [] newElement;
3572	delete [] mark;
3573	// get feasible
3574	if (whichWay < `0.0`) {
3575	newModel.setOptimizationDirection(`1.0`);
3576	double * objective = newModel.objective();
3577	for (int iColumn = `0`; iColumn < numberColumns_; iColumn++)
3578	objective[iColumn] = - objective[iColumn];
3579	}
3580	newModel.primal(`1`);
3581	// still infeasible
3582	if (newModel.problemStatus() == `1`) {
3583	delete [] listNonLinearColumn;
3584	return `0`;
3585	} else if (newModel.problemStatus() == `2`) {
3586	// unbounded - add bounds
3587	double * columnLower = newModel.columnLower();
3588	double * columnUpper = newModel.columnUpper();
3589	for (int i = `0`; i < numberColumns_; i++) {
3590	columnLower[i] = CoinMax(-`1.0e8`, columnLower[i]);
3591	columnUpper[i] = CoinMin(`1.0e8`, columnUpper[i]);
3592	}
3593	newModel.primal(`1`);
3594	}
3595	int numberRows = newModel.numberRows();
3596	double * columnLower = newModel.columnLower();
3597	double * columnUpper = newModel.columnUpper();
3598	double * objective = newModel.objective();
3599	double * rowLower = newModel.rowLower();
3600	double * rowUpper = newModel.rowUpper();
3601	double * solution = newModel.primalColumnSolution();
3602	int jNon;
3603	int * last[`3`];
3604
3605	double * trust = new double[numberNonLinearColumns];
3606	double * trueLower = new double[numberNonLinearColumns];
3607	double * trueUpper = new double[numberNonLinearColumns];
3608	double objectiveOffset;
3609	double objectiveOffset2;
3610	getDblParam(ClpObjOffset, objectiveOffset);
3611	objectiveOffset *= whichWay;
3612	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3613	iColumn = listNonLinearColumn[jNon];
3614	double upper = columnUpper[iColumn];
3615	double lower = columnLower[iColumn];
3616	if (solution[iColumn] < lower)
3617	solution[iColumn] = lower;
3618	else if (solution[iColumn] > upper)
3619	solution[iColumn] = upper;
3620	#if 0
3621	double large = CoinMax(`1000.0`, `10.0` * fabs(solution[iColumn]));
3622	if (upper > `1.0e10`)
3623	upper = solution[iColumn] + large;
3624	if (lower < -`1.0e10`)
3625	lower = solution[iColumn] - large;
3626	#else
3627	upper = solution[iColumn] + `0.5`;
3628	lower = solution[iColumn] - `0.5`;
3629	#endif
3630	//columnUpper[iColumn]=upper;
3631	trust[jNon] = `0.05` * (`1.0` + upper - lower);
3632	trueLower[jNon] = columnLower[iColumn];
3633	//trueUpper[jNon]=upper;
3634	trueUpper[jNon] = columnUpper[iColumn];
3635	}
3636	bool tryFix = false;
3637	int iPass;
3638	double lastObjective = `1.0e31`;
3639	double lastGoodObjective = `1.0e31`;
3640	double * bestSolution = NULL;
3641	double * saveSolution = new double [numberColumns2+numberRows];
3642	char * saveStatus = new char [numberColumns2+numberRows];
3643	double targetDrop = `1.0e31`;
3644	// 1 bound up, 2 up, -1 bound down, -2 down, 0 no change
3645	for (iPass = `0`; iPass < `3`; iPass++) {
3646	last[iPass] = new int[numberNonLinearColumns];
3647	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++)
3648	last[iPass][jNon] = `0`;
3649	}
3650	int numberZeroPasses = `0`;
3651	bool zeroTargetDrop = false;
3652	double * gradient = new double [numberColumns_];
3653	bool goneFeasible = false;
3654	// keep sum of artificials
3655	#define KEEP_SUM 5
3656	double sumArt[KEEP_SUM];
3657	for (jNon = `0`; jNon < KEEP_SUM; jNon++)
3658	sumArt[jNon] = COIN_DBL_MAX;
3659	#define SMALL_FIX 0.0
3660	for (iPass = `0`; iPass < numberPasses; iPass++) {
3661	objectiveOffset2 = objectiveOffset;
3662	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3663	iColumn = listNonLinearColumn[jNon];
3664	if (solution[iColumn] < trueLower[jNon])
3665	solution[iColumn] = trueLower[jNon];
3666	else if (solution[iColumn] > trueUpper[jNon])
3667	solution[iColumn] = trueUpper[jNon];
3668	columnLower[iColumn] = CoinMax(solution[iColumn]
3669	- trust[jNon],
3670	trueLower[jNon]);
3671	if (!trueLower[jNon] && columnLower[iColumn] < SMALL_FIX)
3672	columnLower[iColumn] = SMALL_FIX;
3673	columnUpper[iColumn] = CoinMin(solution[iColumn]
3674	+ trust[jNon],
3675	trueUpper[jNon]);
3676	if (!trueLower[jNon])
3677	columnUpper[iColumn] = CoinMax(columnUpper[iColumn],
3678	columnLower[iColumn] + SMALL_FIX);
3679	if (!trueLower[jNon] && tryFix &&
3680	columnLower[iColumn] == SMALL_FIX &&
3681	columnUpper[iColumn] < `3.0` * SMALL_FIX) {
3682	columnLower[iColumn] = `0.0`;
3683	solution[iColumn] = `0.0`;
3684	columnUpper[iColumn] = `0.0`;
3685	printf("fixing %d\n", iColumn);
3686	}
3687	}
3688	// redo matrix
3689	double offset;
3690	CoinPackedMatrix newMatrix(saveMatrix);
3691	// get matrix data pointers
3692	column = newMatrix.getIndices();
3693	rowStart = newMatrix.getVectorStarts();
3694	rowLength = newMatrix.getVectorLengths();
3695	// make sure x updated
3696	if (numberConstraints)
3697	constraints[`0`]->newXValues();
3698	else
3699	trueObjective->newXValues();
3700	double * changeableElement = newMatrix.getMutableElements();
3701	if (trueObjective) {
3702	CoinMemcpyN(trueObjective->gradient(this, solution, offset, true, `2`), numberColumns_,
3703	objective);
3704	} else {
3705	CoinMemcpyN(objective_->gradient(this, solution, offset, true, `2`), numberColumns_,
3706	objective);
3707	}
3708	if (whichWay < `0.0`) {
3709	for (int iColumn = `0`; iColumn < numberColumns_; iColumn++)
3710	objective[iColumn] = - objective[iColumn];
3711	}
3712	for (iConstraint = `0`; iConstraint < numberConstraints; iConstraint++) {
3713	ClpConstraint * constraint = constraints[iConstraint];
3714	int iRow = constraint->rowNumber();
3715	double functionValue;
3716	#ifndef NDEBUG
3717	int numberErrors =
3718	#endif
3719	constraint->gradient(&newModel, solution, gradient, functionValue, offset);
3720	assert (!numberErrors);
3721	// double dualValue = newModel.dualRowSolution()[iRow];
3722	int numberCoefficients = constraint->numberCoefficients();
3723	for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + numberCoefficients; j++) {
3724	int iColumn = column[j];
3725	changeableElement[j] = gradient[iColumn];
3726	//objective[iColumn] -= dualValuegradient[iColumn];*
3727	gradient[iColumn] = `0.0`;
3728	}
3729	for (int k = `0`; k < numberColumns_; k++)
3730	assert (!gradient[k]);
3731	if (rowLower_[iRow] > -`1.0e20`)
3732	rowLower[iRow] = rowLower_[iRow] - offset;
3733	if (rowUpper_[iRow] < `1.0e20`)
3734	rowUpper[iRow] = rowUpper_[iRow] - offset;
3735	}
3736	// Replace matrix
3737	// Get a column copy in standard format
3738	CoinPackedMatrix * columnCopy = new CoinPackedMatrix ();
3739	columnCopy->reverseOrderedCopyOf(newMatrix);
3740	newModel.replaceMatrix(columnCopy, true);
3741	// solve
3742	newModel.primal(`1`);
3743	if (newModel.status() == `1`) {
3744	// Infeasible!
3745	newModel.allSlackBasis();
3746	newModel.primal();
3747	newModel.writeMps("infeas.mps");
3748	assert(!newModel.status());
3749	}
3750	double sumInfeas = `0.0`;
3751	int numberInfeas = `0`;
3752	for (iColumn = numberColumns_; iColumn < numberColumns2; iColumn++) {
3753	if (solution[iColumn] > `1.0e-8`) {
3754	numberInfeas++;
3755	sumInfeas += solution[iColumn];
3756	}
3757	}
3758	printf("%d artificial infeasibilities - summing to %g\n",
3759	numberInfeas, sumInfeas);
3760	for (jNon = `0`; jNon < KEEP_SUM - `1`; jNon++)
3761	sumArt[jNon] = sumArt[jNon+`1`];
3762	sumArt[KEEP_SUM-`1`] = sumInfeas;
3763	if (sumInfeas > `0.01` && sumInfeas * `1.1` > sumArt[`0`] && penalties[`1`] < `1.0e7`) {
3764	// not doing very well - increase - be more sophisticated later
3765	lastObjective = COIN_DBL_MAX;
3766	for (jNon = `0`; jNon < KEEP_SUM; jNon++)
3767	sumArt[jNon] = COIN_DBL_MAX;
3768	//for (iColumn=numberColumns_;iColumn<numberColumns2;iColumn+=SEGMENTS) {
3769	//objective[iColumn+1] = 1.5;*
3770	//}
3771	penalties[`1`] *= `1.5`;
3772	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++)
3773	if (trust[jNon] > `0.1`)
3774	trust[jNon] *= `2.0`;
3775	else
3776	trust[jNon] = `0.1`;
3777	}
3778	if (sumInfeas < `0.001` && !goneFeasible) {
3779	goneFeasible = true;
3780	penalties[`0`] = `1.0e-3`;
3781	penalties[`1`] = `1.0e6`;
3782	printf("Got feasible\n");
3783	}
3784	double infValue = `0.0`;
3785	double objValue = `0.0`;
3786	// make sure x updated
3787	if (numberConstraints)
3788	constraints[`0`]->newXValues();
3789	else
3790	trueObjective->newXValues();
3791	if (trueObjective) {
3792	objValue = trueObjective->objectiveValue(this, solution);
3793	printf("objective offset %g\n", offset);
3794	objectiveOffset2 = objectiveOffset + offset; // ? sign
3795	newModel.setObjectiveOffset(objectiveOffset2);
3796	} else {
3797	objValue = objective_->objectiveValue(this, solution);
3798	}
3799	objValue *= whichWay;
3800	double infPenalty = `0.0`;
3801	// This penalty is for target drop
3802	double infPenalty2 = `0.0`;
3803	//const int row = columnCopy->getIndices();*
3804	//const CoinBigIndex columnStart = columnCopy->getVectorStarts();*
3805	//const int columnLength = columnCopy->getVectorLengths();*
3806	//const double element = columnCopy->getElements();*
3807	double * cost = newModel.objective();
3808	column = newMatrix.getIndices();
3809	rowStart = newMatrix.getVectorStarts();
3810	rowLength = newMatrix.getVectorLengths();
3811	elementByRow = newMatrix.getElements();
3812	int jColumn = numberColumns_;
3813	double objectiveAdjustment = `0.0`;
3814	for (iConstraint = `0`; iConstraint < numberConstraints; iConstraint++) {
3815	ClpConstraint * constraint = constraints[iConstraint];
3816	int iRow = constraint->rowNumber();
3817	double functionValue = constraint->functionValue(this, solution);
3818	double dualValue = newModel.dualRowSolution()[iRow];
3819	if (numberConstraints < -`50`)
3820	printf("For row %d current value is %g (row activity %g) , dual is %g\n", iRow, functionValue,
3821	newModel.primalRowSolution()[iRow],
3822	dualValue);
3823	double movement = newModel.primalRowSolution()[iRow] + constraint->offset();
3824	movement = fabs((movement - functionValue) * dualValue);
3825	infPenalty2 += movement;
3826	double sumOfActivities = `0.0`;
3827	for (CoinBigIndex j = rowStart[iRow]; j < rowStart[iRow] + rowLength[iRow]; j++) {
3828	int iColumn = column[j];
3829	sumOfActivities += fabs(solution[iColumn] * elementByRow[j]);
3830	}
3831	if (rowLower_[iRow] > -`1.0e20`) {
3832	if (functionValue < rowLower_[iRow] - `1.0e-5`) {
3833	double infeasibility = rowLower_[iRow] - functionValue;
3834	double thisPenalty = `0.0`;
3835	infValue += infeasibility;
3836	double boundMultiplier = `1.0`;
3837	if (sumOfActivities < `0.001`)
3838	boundMultiplier = `0.1`;
3839	else if (sumOfActivities > `100.0`)
3840	boundMultiplier = `10.0`;
3841	int k;
3842	assert (dualValue >= -`1.0e-5`);
3843	dualValue = CoinMax(dualValue, `0.0`);
3844	for ( k = `0`; k < SEGMENTS; k++) {
3845	if (infeasibility <= `0`)
3846	break;
3847	double thisPart = CoinMin(infeasibility, bounds[k]);
3848	thisPenalty += thisPart * cost[jColumn+k];
3849	infeasibility -= thisPart;
3850	}
3851	infeasibility = functionValue - rowUpper_[iRow];
3852	double newPenalty = `0.0`;
3853	for ( k = `0`; k < SEGMENTS; k++) {
3854	double thisPart = CoinMin(infeasibility, bounds[k]);
3855	cost[jColumn+k] = CoinMax(penalties[k], dualValue + `1.0e-3`);
3856	newPenalty += thisPart * cost[jColumn+k];
3857	infeasibility -= thisPart;
3858	}
3859	infPenalty += thisPenalty;
3860	objectiveAdjustment += CoinMax(`0.0`, newPenalty - thisPenalty);
3861	}
3862	jColumn += SEGMENTS;
3863	}
3864	if (rowUpper_[iRow] < `1.0e20`) {
3865	if (functionValue > rowUpper_[iRow] + `1.0e-5`) {
3866	double infeasibility = functionValue - rowUpper_[iRow];
3867	double thisPenalty = `0.0`;
3868	infValue += infeasibility;
3869	double boundMultiplier = `1.0`;
3870	if (sumOfActivities < `0.001`)
3871	boundMultiplier = `0.1`;
3872	else if (sumOfActivities > `100.0`)
3873	boundMultiplier = `10.0`;
3874	int k;
3875	dualValue = -dualValue;
3876	assert (dualValue >= -`1.0e-5`);
3877	dualValue = CoinMax(dualValue, `0.0`);
3878	for ( k = `0`; k < SEGMENTS; k++) {
3879	if (infeasibility <= `0`)
3880	break;
3881	double thisPart = CoinMin(infeasibility, bounds[k]);
3882	thisPenalty += thisPart * cost[jColumn+k];
3883	infeasibility -= thisPart;
3884	}
3885	infeasibility = functionValue - rowUpper_[iRow];
3886	double newPenalty = `0.0`;
3887	for ( k = `0`; k < SEGMENTS; k++) {
3888	double thisPart = CoinMin(infeasibility, bounds[k]);
3889	cost[jColumn+k] = CoinMax(penalties[k], dualValue + `1.0e-3`);
3890	newPenalty += thisPart * cost[jColumn+k];
3891	infeasibility -= thisPart;
3892	}
3893	infPenalty += thisPenalty;
3894	objectiveAdjustment += CoinMax(`0.0`, newPenalty - thisPenalty);
3895	}
3896	jColumn += SEGMENTS;
3897	}
3898	}
3899	// adjust last objective value
3900	lastObjective += objectiveAdjustment;
3901	if (infValue)
3902	printf("Sum infeasibilities %g - penalty %g ", infValue, infPenalty);
3903	if (objectiveOffset2)
3904	printf("offset2 %g ", objectiveOffset2);
3905	objValue -= objectiveOffset2;
3906	printf("True objective %g or maybe %g (with penalty %g) -pen2 %g %g\n", objValue,
3907	objValue + objectiveOffset2, objValue + objectiveOffset2 + infPenalty, infPenalty2, penalties[`1`]);
3908	double useObjValue = objValue + objectiveOffset2 + infPenalty;
3909	objValue += infPenalty + infPenalty2;
3910	objValue = useObjValue;
3911	if (iPass) {
3912	double drop = lastObjective - objValue;
3913	std::cout << "True drop was " << drop << std::endl;
3914	if (drop < -`0.05` * fabs(objValue) - `1.0e-4`) {
3915	// pretty bad - go back and halve
3916	CoinMemcpyN(saveSolution, numberColumns2, solution);
3917	CoinMemcpyN(saveSolution + numberColumns2,
3918	numberRows, newModel.primalRowSolution());
3919	CoinMemcpyN(reinterpret_cast<unsigned char *> (saveStatus),
3920	numberColumns2 + numberRows, newModel.statusArray());
3921	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++)
3922	if (trust[jNon] > `0.1`)
3923	trust[jNon] *= `0.5`;
3924	else
3925	trust[jNon] *= `0.9`;
3926
3927	printf("** Halving trust\n");
3928	objValue = lastObjective;
3929	continue;
3930	} else if ((iPass % `25`) == -`1`) {
3931	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++)
3932	trust[jNon] *= `2.0`;
3933	printf("** Doubling trust\n");
3934	}
3935	int * temp = last[`2`];
3936	last[`2`] = last[`1`];
3937	last[`1`] = last[`0`];
3938	last[`0`] = temp;
3939	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3940	iColumn = listNonLinearColumn[jNon];
3941	double change = solution[iColumn] - saveSolution[iColumn];
3942	if (change < -`1.0e-5`) {
3943	if (fabs(change + trust[jNon]) < `1.0e-5`)
3944	temp[jNon] = -`1`;
3945	else
3946	temp[jNon] = -`2`;
3947	} else if(change > `1.0e-5`) {
3948	if (fabs(change - trust[jNon]) < `1.0e-5`)
3949	temp[jNon] = `1`;
3950	else
3951	temp[jNon] = `2`;
3952	} else {
3953	temp[jNon] = `0`;
3954	}
3955	}
3956	double maxDelta = `0.0`;
3957	double smallestTrust = `1.0e31`;
3958	double smallestNonLinearGap = `1.0e31`;
3959	double smallestGap = `1.0e31`;
3960	bool increasing = false;
3961	for (iColumn = `0`; iColumn < numberColumns_; iColumn++) {
3962	double gap = columnUpper[iColumn] - columnLower[iColumn];
3963	assert (gap >= `0.0`);
3964	if (gap)
3965	smallestGap = CoinMin(smallestGap, gap);
3966	}
3967	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
3968	iColumn = listNonLinearColumn[jNon];
3969	double gap = columnUpper[iColumn] - columnLower[iColumn];
3970	assert (gap >= `0.0`);
3971	if (gap) {
3972	smallestNonLinearGap = CoinMin(smallestNonLinearGap, gap);
3973	if (gap < `1.0e-7` && iPass == `1`) {
3974	printf("Small gap %d %d %g %g %g\n",
3975	jNon, iColumn, columnLower[iColumn], columnUpper[iColumn],
3976	gap);
3977	//trueUpper[jNon]=trueLower[jNon];
3978	//columnUpper[iColumn]=columnLower[iColumn];
3979	}
3980	}
3981	maxDelta = CoinMax(maxDelta,
3982	fabs(solution[iColumn] - saveSolution[iColumn]));
3983	if (last[`0`][jNon]*last[`1`][jNon] < `0`) {
3984	// halve
3985	if (trust[jNon] > `1.0`)
3986	trust[jNon] *= `0.5`;
3987	else
3988	trust[jNon] *= `0.7`;
3989	} else {
3990	// ? only increase if +=1 ?
3991	if (last[`0`][jNon] == last[`1`][jNon] &&
3992	last[`0`][jNon] == last[`2`][jNon] &&
3993	last[`0`][jNon]) {
3994	trust[jNon] *= `1.8`;
3995	increasing = true;
3996	}
3997	}
3998	smallestTrust = CoinMin(smallestTrust, trust[jNon]);
3999	}
4000	std::cout << "largest delta is " << maxDelta
4001	<< ", smallest trust is " << smallestTrust
4002	<< ", smallest gap is " << smallestGap
4003	<< ", smallest nonlinear gap is " << smallestNonLinearGap
4004	<< std::endl;
4005	if (iPass > `200`) {
4006	//double useObjValue = objValue+objectiveOffset2+infPenalty;
4007	if (useObjValue + `1.0e-4` > lastGoodObjective && iPass > `250`) {
4008	std::cout << "Exiting as objective not changing much" << std::endl;
4009	break;
4010	} else if (useObjValue < lastGoodObjective) {
4011	lastGoodObjective = useObjValue;
4012	if (!bestSolution)
4013	bestSolution = new double [numberColumns2];
4014	CoinMemcpyN(solution, numberColumns2, bestSolution);
4015	}
4016	}
4017	if (maxDelta < deltaTolerance && !increasing && iPass > `100`) {
4018	numberZeroPasses++;
4019	if (numberZeroPasses == `3`) {
4020	if (tryFix) {
4021	std::cout << "Exiting" << std::endl;
4022	break;
4023	} else {
4024	tryFix = true;
4025	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++)
4026	trust[jNon] = CoinMax(trust[jNon], `1.0e-1`);
4027	numberZeroPasses = `0`;
4028	}
4029	}
4030	} else {
4031	numberZeroPasses = `0`;
4032	}
4033	}
4034	CoinMemcpyN(solution, numberColumns2, saveSolution);
4035	CoinMemcpyN(newModel.primalRowSolution(),
4036	numberRows, saveSolution + numberColumns2);
4037	CoinMemcpyN(newModel.statusArray(),
4038	numberColumns2 + numberRows,
4039	reinterpret_cast<unsigned char *> (saveStatus));
4040
4041	targetDrop = infPenalty + infPenalty2;
4042	if (iPass) {
4043	// get reduced costs
4044	const double * pi = newModel.dualRowSolution();
4045	newModel.matrix()->transposeTimes(pi,
4046	newModel.dualColumnSolution());
4047	double * r = newModel.dualColumnSolution();
4048	for (iColumn = `0`; iColumn < numberColumns_; iColumn++)
4049	r[iColumn] = objective[iColumn] - r[iColumn];
4050	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
4051	iColumn = listNonLinearColumn[jNon];
4052	double dj = r[iColumn];
4053	if (dj < -`1.0e-6`) {
4054	double drop = -dj * (columnUpper[iColumn] - solution[iColumn]);
4055	//double upper = CoinMin(trueUpper[jNon],solution[iColumn]+0.1);
4056	//double drop2 = -dj(upper-solution[iColumn]);*
4057	#if 0
4058	if (drop > `1.0e8` \|\| drop2 > `100.0` * drop \|\| (drop > `1.0e-2` && iPass > `100`))
4059	printf("Big drop %d %g %g %g %g T %g\n",
4060	iColumn, columnLower[iColumn], solution[iColumn],
4061	columnUpper[iColumn], dj, trueUpper[jNon]);
4062	#endif
4063	targetDrop += drop;
4064	if (dj < -`1.0e-1` && trust[jNon] < `1.0e-3`
4065	&& trueUpper[jNon] - solution[iColumn] > `1.0e-2`) {
4066	trust[jNon] *= `1.5`;
4067	//printf("Increasing trust on %d to %g\n",
4068	// iColumn,trust[jNon]);
4069	}
4070	} else if (dj > `1.0e-6`) {
4071	double drop = -dj * (columnLower[iColumn] - solution[iColumn]);
4072	//double lower = CoinMax(trueLower[jNon],solution[iColumn]-0.1);
4073	//double drop2 = -dj(lower-solution[iColumn]);*
4074	#if 0
4075	if (drop > `1.0e8` \|\| drop2 > `100.0` * drop \|\| (drop > `1.0e-2`))
4076	printf("Big drop %d %g %g %g %g T %g\n",
4077	iColumn, columnLower[iColumn], solution[iColumn],
4078	columnUpper[iColumn], dj, trueLower[jNon]);
4079	#endif
4080	targetDrop += drop;
4081	if (dj > `1.0e-1` && trust[jNon] < `1.0e-3`
4082	&& solution[iColumn] - trueLower[jNon] > `1.0e-2`) {
4083	trust[jNon] *= `1.5`;
4084	printf("Increasing trust on %d to %g\n",
4085	iColumn, trust[jNon]);
4086	}
4087	}
4088	}
4089	}
4090	std::cout << "Pass - " << iPass
4091	<< ", target drop is " << targetDrop
4092	<< std::endl;
4093	if (iPass > `1` && targetDrop < `1.0e-5` && zeroTargetDrop)
4094	break;
4095	if (iPass > `1` && targetDrop < `1.0e-5`)
4096	zeroTargetDrop = true;
4097	else
4098	zeroTargetDrop = false;
4099	//if (iPass==5)
4100	//newModel.setLogLevel(63);
4101	lastObjective = objValue;
4102	// take out when ClpPackedMatrix changed
4103	//newModel.scaling(false);
4104	#if 0
4105	CoinMpsIO writer;
4106	writer.setMpsData(*newModel.matrix(), COIN_DBL_MAX,
4107	newModel.getColLower(), newModel.getColUpper(),
4108	newModel.getObjCoefficients(),
4109	(const char) `0` /integrality/*,
4110	newModel.getRowLower(), newModel.getRowUpper(),
4111	NULL, NULL);
4112	writer.writeMps("xx.mps");
4113	#endif
4114	}
4115	delete [] saveSolution;
4116	delete [] saveStatus;
4117	for (iPass = `0`; iPass < `3`; iPass++)
4118	delete [] last[iPass];
4119	for (jNon = `0`; jNon < numberNonLinearColumns; jNon++) {
4120	iColumn = listNonLinearColumn[jNon];
4121	columnLower[iColumn] = trueLower[jNon];
4122	columnUpper[iColumn] = trueUpper[jNon];
4123	}
4124	delete [] trust;
4125	delete [] trueUpper;
4126	delete [] trueLower;
4127	objectiveValue_ = newModel.objectiveValue();
4128	if (bestSolution) {
4129	CoinMemcpyN(bestSolution, numberColumns2, solution);
4130	delete [] bestSolution;
4131	printf("restoring objective of %g\n", lastGoodObjective);
4132	objectiveValue_ = lastGoodObjective;
4133	}
4134	// Simplest way to get true row activity ?
4135	double * rowActivity = newModel.primalRowSolution();
4136	for (iRow = `0`; iRow < numberRows; iRow++) {
4137	double difference;
4138	if (fabs(rowLower_[iRow]) < fabs(rowUpper_[iRow]))
4139	difference = rowLower_[iRow] - rowLower[iRow];
4140	else
4141	difference = rowUpper_[iRow] - rowUpper[iRow];
4142	rowLower[iRow] = rowLower_[iRow];
4143	rowUpper[iRow] = rowUpper_[iRow];
4144	if (difference) {
4145	if (numberRows < `50`)
4146	printf("For row %d activity changes from %g to %g\n",
4147	iRow, rowActivity[iRow], rowActivity[iRow] + difference);
4148	rowActivity[iRow] += difference;
4149	}
4150	}
4151	delete [] listNonLinearColumn;
4152	delete [] gradient;
4153	printf("solution still in newModel - do objective etc!\n");
4154	numberIterations_ = newModel.numberIterations();
4155	problemStatus_ = newModel.problemStatus();
4156	secondaryStatus_ = newModel.secondaryStatus();
4157	CoinMemcpyN(newModel.primalColumnSolution(), numberColumns_, columnActivity_);
4158	// should do status region
4159	CoinZeroN(rowActivity_, numberRows_);
4160	matrix_->times(`1.0`, columnActivity_, rowActivity_);
4161	return `0`;
4162	}
4163

Browse the source code of OGDF/src/coin/Clp/ClpSimplexNonlinear.cpp