Context Navigation

source: trunk/Clp/src/ClpSimplexNonlinear.cpp @ 2385

Last change on this file since 2385 was 2385, checked in by unxusr, 4 months ago
formatting
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 152.2 KB

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

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

Download in other formats: