Context Navigation

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

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago
Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.
Property svn:eol-style set to `native` Property svn:keywords set to `Id`
File size: 184.6 KB

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

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

Download in other formats: