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

1	/ $Id: ClpPackedMatrix.cpp 1864 2012-06-28 10:27:20Z forrest $ /
2	// Copyright (C) 2002, International Business Machines
3	// Corporation and others. All Rights Reserved.
4	// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7
8	#include <cstdio>
9
10	#include "CoinPragma.hpp"
11	#include "CoinIndexedVector.hpp"
12	#include "CoinHelperFunctions.hpp"
13	//#define THREAD
14
15	#include "ClpSimplex.hpp"
16	#include "ClpSimplexDual.hpp"
17	#include "ClpFactorization.hpp"
18	#ifndef SLIM_CLP
19	#include "ClpQuadraticObjective.hpp"
20	#endif
21	// at end to get min/max!
22	#include "ClpPackedMatrix.hpp"
23	#include "ClpMessage.hpp"
24	#ifdef INTEL_MKL
25	#include "mkl_spblas.h"
26	#endif
27
28	//=============================================================================
29	#ifdef COIN_PREFETCH
30	#if 1
31	#define coin_prefetch(mem) \
32	__asm__ __volatile__ ("prefetchnta %0" : : "m" ((reinterpret_cast<char >(mem))))
33	#define coin_prefetch_const(mem) \
34	__asm__ __volatile__ ("prefetchnta %0" : : "m" ((reinterpret_cast<const char >(mem))))
35	#else
36	#define coin_prefetch(mem) \
37	__asm__ __volatile__ ("prefetch %0" : : "m" ((reinterpret_cast<char >(mem))))
38	#define coin_prefetch_const(mem) \
39	__asm__ __volatile__ ("prefetch %0" : : "m" ((reinterpret_cast<const char >(mem))))
40	#endif
41	#else
42	// dummy
43	#define coin_prefetch(mem)
44	#define coin_prefetch_const(mem)
45	#endif
46
47	//#############################################################################
48	// Constructors / Destructor / Assignment
49	//#############################################################################
50
51	//-------------------------------------------------------------------
52	// Default Constructor
53	//-------------------------------------------------------------------
54	ClpPackedMatrix::ClpPackedMatrix ()
55	: ClpMatrixBase (),
56	matrix_(NULL),
57	numberActiveColumns_(`0`),
58	flags_(`2`),
59	rowCopy_(NULL),
60	columnCopy_(NULL)
61	{
62	setType(`1`);
63	}
64
65	//-------------------------------------------------------------------
66	// Copy constructor
67	//-------------------------------------------------------------------
68	ClpPackedMatrix::ClpPackedMatrix (const ClpPackedMatrix & rhs)
69	: ClpMatrixBase (rhs)
70	{
71	#ifndef COIN_SPARSE_MATRIX
72	matrix_ = new CoinPackedMatrix (*(rhs.matrix_), -`1`, -`1`);
73	#else
74	matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -`0`, -`0`);
75	#endif
76	numberActiveColumns_ = rhs.numberActiveColumns_;
77	flags_ = rhs.flags_ & (~`2`);
78	int numberRows = matrix_->getNumRows();
79	if (rhs.rhsOffset_ && numberRows) {
80	rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
81	} else {
82	rhsOffset_ = NULL;
83	}
84	if (rhs.rowCopy_) {
85	assert ((flags_ & `4`) != `0`);
86	rowCopy_ = new ClpPackedMatrix2 (*rhs.rowCopy_);
87	} else {
88	rowCopy_ = NULL;
89	}
90	if (rhs.columnCopy_) {
91	assert ((flags_&(`8` + `16`)) == `8` + `16`);
92	columnCopy_ = new ClpPackedMatrix3 (*rhs.columnCopy_);
93	} else {
94	columnCopy_ = NULL;
95	}
96	}
97	//-------------------------------------------------------------------
98	// assign matrix (for space reasons)
99	//-------------------------------------------------------------------
100	ClpPackedMatrix::ClpPackedMatrix (CoinPackedMatrix * rhs)
101	: ClpMatrixBase ()
102	{
103	matrix_ = rhs;
104	flags_ = matrix_->hasGaps() ? `2` : `0`;
105	numberActiveColumns_ = matrix_->getNumCols();
106	rowCopy_ = NULL;
107	columnCopy_ = NULL;
108	setType(`1`);
109
110	}
111
112	ClpPackedMatrix::ClpPackedMatrix (const CoinPackedMatrix & rhs)
113	: ClpMatrixBase ()
114	{
115	#ifndef COIN_SPARSE_MATRIX
116	matrix_ = new CoinPackedMatrix (rhs, -`1`, -`1`);
117	#else
118	matrix_ = new CoinPackedMatrix(rhs, -`0`, -`0`);
119	#endif
120	numberActiveColumns_ = matrix_->getNumCols();
121	rowCopy_ = NULL;
122	flags_ = `0`;
123	columnCopy_ = NULL;
124	setType(`1`);
125
126	}
127
128	//-------------------------------------------------------------------
129	// Destructor
130	//-------------------------------------------------------------------
131	ClpPackedMatrix::~ClpPackedMatrix ()
132	{
133	delete matrix_;
134	delete rowCopy_;
135	delete columnCopy_;
136	}
137
138	//----------------------------------------------------------------
139	// Assignment operator
140	//-------------------------------------------------------------------
141	ClpPackedMatrix &
142	ClpPackedMatrix::operator=(const ClpPackedMatrix& rhs)
143	{
144	if (this != &rhs) {
145	ClpMatrixBase::operator=(rhs);
146	delete matrix_;
147	#ifndef COIN_SPARSE_MATRIX
148	matrix_ = new CoinPackedMatrix (*(rhs.matrix_));
149	#else
150	matrix_ = new CoinPackedMatrix(*(rhs.matrix_), -`0`, -`0`);
151	#endif
152	numberActiveColumns_ = rhs.numberActiveColumns_;
153	flags_ = rhs.flags_;
154	delete rowCopy_;
155	delete columnCopy_;
156	if (rhs.rowCopy_) {
157	assert ((flags_ & `4`) != `0`);
158	rowCopy_ = new ClpPackedMatrix2 (*rhs.rowCopy_);
159	} else {
160	rowCopy_ = NULL;
161	}
162	if (rhs.columnCopy_) {
163	assert ((flags_&(`8` + `16`)) == `8` + `16`);
164	columnCopy_ = new ClpPackedMatrix3 (*rhs.columnCopy_);
165	} else {
166	columnCopy_ = NULL;
167	}
168	}
169	return *this;
170	}
171	//-------------------------------------------------------------------
172	// Clone
173	//-------------------------------------------------------------------
174	ClpMatrixBase * ClpPackedMatrix::clone() const
175	{
176	return new ClpPackedMatrix (*this);
177	}
178	// Copy contents - resizing if necessary - otherwise re-use memory
179	void
180	ClpPackedMatrix::copy(const ClpPackedMatrix * rhs)
181	{
182	//this = rhs;
183	assert(numberActiveColumns_ == rhs->numberActiveColumns_);
184	assert (matrix_->isColOrdered() == rhs->matrix_->isColOrdered());
185	matrix_->copyReuseArrays(*rhs->matrix_);
186	}
187	/ Subset clone (without gaps). Duplicates are allowed*
188	and order is as given /*
189	ClpMatrixBase *
190	ClpPackedMatrix::subsetClone (int numberRows, const int * whichRows,
191	int numberColumns,
192	const int * whichColumns) const
193	{
194	return new ClpPackedMatrix (*this, numberRows, whichRows,
195	numberColumns, whichColumns);
196	}
197	/ Subset constructor (without gaps). Duplicates are allowed*
198	and order is as given /*
199	ClpPackedMatrix::ClpPackedMatrix (
200	const ClpPackedMatrix & rhs,
201	int numberRows, const int * whichRows,
202	int numberColumns, const int * whichColumns)
203	: ClpMatrixBase (rhs)
204	{
205	matrix_ = new CoinPackedMatrix (*(rhs.matrix_), numberRows, whichRows,
206	numberColumns, whichColumns);
207	numberActiveColumns_ = matrix_->getNumCols();
208	rowCopy_ = NULL;
209	flags_ = rhs.flags_ & (~`2`); // no gaps
210	columnCopy_ = NULL;
211	}
212	ClpPackedMatrix::ClpPackedMatrix (
213	const CoinPackedMatrix & rhs,
214	int numberRows, const int * whichRows,
215	int numberColumns, const int * whichColumns)
216	: ClpMatrixBase ()
217	{
218	matrix_ = new CoinPackedMatrix (rhs, numberRows, whichRows,
219	numberColumns, whichColumns);
220	numberActiveColumns_ = matrix_->getNumCols();
221	rowCopy_ = NULL;
222	flags_ = `2`;
223	columnCopy_ = NULL;
224	setType(`1`);
225	}
226
227	/ Returns a new matrix in reverse order without gaps /
228	ClpMatrixBase *
229	ClpPackedMatrix::reverseOrderedCopy() const
230	{
231	ClpPackedMatrix * copy = new ClpPackedMatrix ();
232	copy->matrix_ = new CoinPackedMatrix ();
233	copy->matrix_->setExtraGap(`0.0`);
234	copy->matrix_->setExtraMajor(`0.0`);
235	copy->matrix_->reverseOrderedCopyOf(*matrix_);
236	//copy->matrix_->removeGaps();
237	copy->numberActiveColumns_ = copy->matrix_->getNumCols();
238	copy->flags_ = flags_ & (~`2`); // no gaps
239	return copy;
240	}
241	//unscaled versions
242	void
243	ClpPackedMatrix::times(double scalar,
244	const double * x, double * y) const
245	{
246	int iRow, iColumn;
247	// get matrix data pointers
248	const int * row = matrix_->getIndices();
249	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
250	const double * elementByColumn = matrix_->getElements();
251	//memset(y,0,matrix_->getNumRows()sizeof(double));*
252	assert (((flags_ & `2`) != `0`) == matrix_->hasGaps());
253	if (!(flags_ & `2`)) {
254	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
255	CoinBigIndex j;
256	double value = x[iColumn];
257	if (value) {
258	CoinBigIndex start = columnStart[iColumn];
259	CoinBigIndex end = columnStart[iColumn+`1`];
260	value *= scalar;
261	for (j = start; j < end; j++) {
262	iRow = row[j];
263	y[iRow] += value * elementByColumn[j];
264	}
265	}
266	}
267	} else {
268	const int * columnLength = matrix_->getVectorLengths();
269	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
270	CoinBigIndex j;
271	double value = x[iColumn];
272	if (value) {
273	CoinBigIndex start = columnStart[iColumn];
274	CoinBigIndex end = start + columnLength[iColumn];
275	value *= scalar;
276	for (j = start; j < end; j++) {
277	iRow = row[j];
278	y[iRow] += value * elementByColumn[j];
279	}
280	}
281	}
282	}
283	}
284	void
285	ClpPackedMatrix::transposeTimes(double scalar,
286	const double * x, double * y) const
287	{
288	int iColumn;
289	// get matrix data pointers
290	const int * row = matrix_->getIndices();
291	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
292	const double * elementByColumn = matrix_->getElements();
293	if (!(flags_ & `2`)) {
294	if (scalar == -`1.0`) {
295	CoinBigIndex start = columnStart[`0`];
296	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
297	CoinBigIndex j;
298	CoinBigIndex next = columnStart[iColumn+`1`];
299	double value = y[iColumn];
300	for (j = start; j < next; j++) {
301	int jRow = row[j];
302	value -= x[jRow] * elementByColumn[j];
303	}
304	start = next;
305	y[iColumn] = value;
306	}
307	} else {
308	CoinBigIndex start = columnStart[`0`];
309	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
310	CoinBigIndex j;
311	CoinBigIndex next = columnStart[iColumn+`1`];
312	double value = `0.0`;
313	for (j = start; j < next; j++) {
314	int jRow = row[j];
315	value += x[jRow] * elementByColumn[j];
316	}
317	start = next;
318	y[iColumn] += value * scalar;
319	}
320	}
321	} else {
322	const int * columnLength = matrix_->getVectorLengths();
323	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
324	CoinBigIndex j;
325	double value = `0.0`;
326	CoinBigIndex start = columnStart[iColumn];
327	CoinBigIndex end = start + columnLength[iColumn];
328	for (j = start; j < end; j++) {
329	int jRow = row[j];
330	value += x[jRow] * elementByColumn[j];
331	}
332	y[iColumn] += value * scalar;
333	}
334	}
335	}
336	void
337	ClpPackedMatrix::times(double scalar,
338	const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
339	const double * COIN_RESTRICT rowScale,
340	const double * COIN_RESTRICT columnScale) const
341	{
342	if (rowScale) {
343	int iRow, iColumn;
344	// get matrix data pointers
345	const int * COIN_RESTRICT row = matrix_->getIndices();
346	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
347	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
348	if (!(flags_ & `2`)) {
349	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
350	double value = x[iColumn];
351	if (value) {
352	// scaled
353	value = scalar columnScale[iColumn];
354	CoinBigIndex start = columnStart[iColumn];
355	CoinBigIndex end = columnStart[iColumn+`1`];
356	CoinBigIndex j;
357	for (j = start; j < end; j++) {
358	iRow = row[j];
359	y[iRow] += value * elementByColumn[j] * rowScale[iRow];
360	}
361	}
362	}
363	} else {
364	const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
365	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
366	double value = x[iColumn];
367	if (value) {
368	// scaled
369	value = scalar columnScale[iColumn];
370	CoinBigIndex start = columnStart[iColumn];
371	CoinBigIndex end = start + columnLength[iColumn];
372	CoinBigIndex j;
373	for (j = start; j < end; j++) {
374	iRow = row[j];
375	y[iRow] += value * elementByColumn[j] * rowScale[iRow];
376	}
377	}
378	}
379	}
380	} else {
381	times(scalar, x, y);
382	}
383	}
384	void
385	ClpPackedMatrix::transposeTimes( double scalar,
386	const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
387	const double * COIN_RESTRICT rowScale,
388	const double * COIN_RESTRICT columnScale,
389	double * COIN_RESTRICT spare) const
390	{
391	if (rowScale) {
392	int iColumn;
393	// get matrix data pointers
394	const int * COIN_RESTRICT row = matrix_->getIndices();
395	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
396	const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
397	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
398	if (!spare) {
399	if (!(flags_ & `2`)) {
400	CoinBigIndex start = columnStart[`0`];
401	if (scalar == -`1.0`) {
402	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
403	CoinBigIndex j;
404	CoinBigIndex next = columnStart[iColumn+`1`];
405	double value = `0.0`;
406	// scaled
407	for (j = start; j < next; j++) {
408	int jRow = row[j];
409	value += x[jRow] * elementByColumn[j] * rowScale[jRow];
410	}
411	start = next;
412	y[iColumn] -= value * columnScale[iColumn];
413	}
414	} else {
415	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
416	CoinBigIndex j;
417	CoinBigIndex next = columnStart[iColumn+`1`];
418	double value = `0.0`;
419	// scaled
420	for (j = start; j < next; j++) {
421	int jRow = row[j];
422	value += x[jRow] * elementByColumn[j] * rowScale[jRow];
423	}
424	start = next;
425	y[iColumn] += value * scalar * columnScale[iColumn];
426	}
427	}
428	} else {
429	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
430	CoinBigIndex j;
431	double value = `0.0`;
432	// scaled
433	for (j = columnStart[iColumn];
434	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
435	int jRow = row[j];
436	value += x[jRow] * elementByColumn[j] * rowScale[jRow];
437	}
438	y[iColumn] += value * scalar * columnScale[iColumn];
439	}
440	}
441	} else {
442	// can use spare region
443	int iRow;
444	int numberRows = matrix_->getNumRows();
445	for (iRow = `0`; iRow < numberRows; iRow++) {
446	double value = x[iRow];
447	if (value)
448	spare[iRow] = value * rowScale[iRow];
449	else
450	spare[iRow] = `0.0`;
451	}
452	if (!(flags_ & `2`)) {
453	CoinBigIndex start = columnStart[`0`];
454	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
455	CoinBigIndex j;
456	CoinBigIndex next = columnStart[iColumn+`1`];
457	double value = `0.0`;
458	// scaled
459	for (j = start; j < next; j++) {
460	int jRow = row[j];
461	value += spare[jRow] * elementByColumn[j];
462	}
463	start = next;
464	y[iColumn] += value * scalar * columnScale[iColumn];
465	}
466	} else {
467	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
468	CoinBigIndex j;
469	double value = `0.0`;
470	// scaled
471	for (j = columnStart[iColumn];
472	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
473	int jRow = row[j];
474	value += spare[jRow] * elementByColumn[j];
475	}
476	y[iColumn] += value * scalar * columnScale[iColumn];
477	}
478	}
479	// no need to zero out
480	//for (iRow=0;iRow<numberRows;iRow++)
481	//spare[iRow] = 0.0;
482	}
483	} else {
484	transposeTimes(scalar, x, y);
485	}
486	}
487	void
488	ClpPackedMatrix::transposeTimesSubset( int number,
489	const int * which,
490	const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
491	const double * COIN_RESTRICT rowScale,
492	const double * COIN_RESTRICT columnScale,
493	double * COIN_RESTRICT spare) const
494	{
495	// get matrix data pointers
496	const int * COIN_RESTRICT row = matrix_->getIndices();
497	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
498	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
499	if (!spare \|\| !rowScale) {
500	if (rowScale) {
501	for (int jColumn = `0`; jColumn < number; jColumn++) {
502	int iColumn = which[jColumn];
503	CoinBigIndex j;
504	CoinBigIndex start = columnStart[iColumn];
505	CoinBigIndex next = columnStart[iColumn+`1`];
506	double value = `0.0`;
507	for (j = start; j < next; j++) {
508	int jRow = row[j];
509	value += x[jRow] * elementByColumn[j] * rowScale[jRow];
510	}
511	y[iColumn] -= value * columnScale[iColumn];
512	}
513	} else {
514	for (int jColumn = `0`; jColumn < number; jColumn++) {
515	int iColumn = which[jColumn];
516	CoinBigIndex j;
517	CoinBigIndex start = columnStart[iColumn];
518	CoinBigIndex next = columnStart[iColumn+`1`];
519	double value = `0.0`;
520	for (j = start; j < next; j++) {
521	int jRow = row[j];
522	value += x[jRow] * elementByColumn[j];
523	}
524	y[iColumn] -= value;
525	}
526	}
527	} else {
528	// can use spare region
529	int iRow;
530	int numberRows = matrix_->getNumRows();
531	for (iRow = `0`; iRow < numberRows; iRow++) {
532	double value = x[iRow];
533	if (value)
534	spare[iRow] = value * rowScale[iRow];
535	else
536	spare[iRow] = `0.0`;
537	}
538	for (int jColumn = `0`; jColumn < number; jColumn++) {
539	int iColumn = which[jColumn];
540	CoinBigIndex j;
541	CoinBigIndex start = columnStart[iColumn];
542	CoinBigIndex next = columnStart[iColumn+`1`];
543	double value = `0.0`;
544	for (j = start; j < next; j++) {
545	int jRow = row[j];
546	value += spare[jRow] * elementByColumn[j];
547	}
548	y[iColumn] -= value * columnScale[iColumn];
549	}
550	}
551	}
552	/ Return <code>x * A + y</code> in <code>z</code>.*
553	Squashes small elements and knows about ClpSimplex /*
554	void
555	ClpPackedMatrix::transposeTimes(const ClpSimplex * model, double scalar,
556	const CoinIndexedVector * rowArray,
557	CoinIndexedVector * y,
558	CoinIndexedVector * columnArray) const
559	{
560	columnArray->clear();
561	double * pi = rowArray->denseVector();
562	int numberNonZero = `0`;
563	int * index = columnArray->getIndices();
564	double * array = columnArray->denseVector();
565	int numberInRowArray = rowArray->getNumElements();
566	// maybe I need one in OsiSimplex
567	double zeroTolerance = model->zeroTolerance();
568	#if 0 //def COIN_DEVELOP
569	if (zeroTolerance != `1.0e-13`) {
570	printf("small element in matrix - zero tolerance %g\n", zeroTolerance);
571	}
572	#endif
573	int numberRows = model->numberRows();
574	ClpPackedMatrix* rowCopy =
575	static_cast< ClpPackedMatrix*>(model->rowCopy());
576	bool packed = rowArray->packedMode();
577	double factor = (numberRows < `100`) ? `0.25` : `0.35`;
578	factor = `0.5`;
579	// We may not want to do by row if there may be cache problems
580	// It would be nice to find L2 cache size - for moment 512K
581	// Be slightly optimistic
582	if (numberActiveColumns_ * sizeof(double) > `1000000`) {
583	if (numberRows * `10` < numberActiveColumns_)
584	factor *= `0.333333333`;
585	else if (numberRows * `4` < numberActiveColumns_)
586	factor *= `0.5`;
587	else if (numberRows * `2` < numberActiveColumns_)
588	factor *= `0.66666666667`;
589	//if (model->numberIterations()%50==0)
590	//printf("%d nonzero\n",numberInRowArray);
591	}
592	// if not packed then bias a bit more towards by column
593	if (!packed)
594	factor *= `0.9`;
595	assert (!y->getNumElements());
596	double multiplierX = `0.8`;
597	double factor2 = factor * multiplierX;
598	if (packed && rowCopy_ && numberInRowArray > `2` && numberInRowArray > factor2 * numberRows &&
599	numberInRowArray < `0.9` * numberRows && `0`) {
600	rowCopy_->transposeTimes(model, rowCopy->matrix_, rowArray, y, columnArray);
601	return;
602	}
603	if (numberInRowArray > factor * numberRows \|\| !rowCopy) {
604	// do by column
605	// If no gaps - can do a bit faster
606	if (!(flags_ & `2`) \|\| columnCopy_) {
607	transposeTimesByColumn( model, scalar,
608	rowArray, y, columnArray);
609	return;
610	}
611	int iColumn;
612	// get matrix data pointers
613	const int * row = matrix_->getIndices();
614	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
615	const int * columnLength = matrix_->getVectorLengths();
616	const double * elementByColumn = matrix_->getElements();
617	const double * rowScale = model->rowScale();
618	#if 0
619	ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
620	if (rowScale && scaledMatrix) {
621	rowScale = NULL;
622	// get matrix data pointers
623	row = scaledMatrix->getIndices();
624	columnStart = scaledMatrix->getVectorStarts();
625	columnLength = scaledMatrix->getVectorLengths();
626	elementByColumn = scaledMatrix->getElements();
627	}
628	#endif
629	if (packed) {
630	// need to expand pi into y
631	assert(y->capacity() >= numberRows);
632	double * piOld = pi;
633	pi = y->denseVector();
634	const int * whichRow = rowArray->getIndices();
635	int i;
636	if (!rowScale) {
637	// modify pi so can collapse to one loop
638	if (scalar == -`1.0`) {
639	for (i = `0`; i < numberInRowArray; i++) {
640	int iRow = whichRow[i];
641	pi[iRow] = -piOld[i];
642	}
643	} else {
644	for (i = `0`; i < numberInRowArray; i++) {
645	int iRow = whichRow[i];
646	pi[iRow] = scalar * piOld[i];
647	}
648	}
649	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
650	double value = `0.0`;
651	CoinBigIndex j;
652	for (j = columnStart[iColumn];
653	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
654	int iRow = row[j];
655	value += pi[iRow] * elementByColumn[j];
656	}
657	if (fabs(value) > zeroTolerance) {
658	array[numberNonZero] = value;
659	index[numberNonZero++] = iColumn;
660	}
661	}
662	} else {
663	#ifdef CLP_INVESTIGATE
664	if (model->clpScaledMatrix())
665	printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
666	#endif
667	// scaled
668	// modify pi so can collapse to one loop
669	if (scalar == -`1.0`) {
670	for (i = `0`; i < numberInRowArray; i++) {
671	int iRow = whichRow[i];
672	pi[iRow] = -piOld[i] * rowScale[iRow];
673	}
674	} else {
675	for (i = `0`; i < numberInRowArray; i++) {
676	int iRow = whichRow[i];
677	pi[iRow] = scalar * piOld[i] * rowScale[iRow];
678	}
679	}
680	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
681	double value = `0.0`;
682	CoinBigIndex j;
683	const double * columnScale = model->columnScale();
684	for (j = columnStart[iColumn];
685	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
686	int iRow = row[j];
687	value += pi[iRow] * elementByColumn[j];
688	}
689	value *= columnScale[iColumn];
690	if (fabs(value) > zeroTolerance) {
691	array[numberNonZero] = value;
692	index[numberNonZero++] = iColumn;
693	}
694	}
695	}
696	// zero out
697	for (i = `0`; i < numberInRowArray; i++) {
698	int iRow = whichRow[i];
699	pi[iRow] = `0.0`;
700	}
701	} else {
702	if (!rowScale) {
703	if (scalar == -`1.0`) {
704	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
705	double value = `0.0`;
706	CoinBigIndex j;
707	for (j = columnStart[iColumn];
708	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
709	int iRow = row[j];
710	value += pi[iRow] * elementByColumn[j];
711	}
712	if (fabs(value) > zeroTolerance) {
713	index[numberNonZero++] = iColumn;
714	array[iColumn] = -value;
715	}
716	}
717	} else {
718	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
719	double value = `0.0`;
720	CoinBigIndex j;
721	for (j = columnStart[iColumn];
722	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
723	int iRow = row[j];
724	value += pi[iRow] * elementByColumn[j];
725	}
726	value *= scalar;
727	if (fabs(value) > zeroTolerance) {
728	index[numberNonZero++] = iColumn;
729	array[iColumn] = value;
730	}
731	}
732	}
733	} else {
734	#ifdef CLP_INVESTIGATE
735	if (model->clpScaledMatrix())
736	printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
737	#endif
738	// scaled
739	if (scalar == -`1.0`) {
740	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
741	double value = `0.0`;
742	CoinBigIndex j;
743	const double * columnScale = model->columnScale();
744	for (j = columnStart[iColumn];
745	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
746	int iRow = row[j];
747	value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
748	}
749	value *= columnScale[iColumn];
750	if (fabs(value) > zeroTolerance) {
751	index[numberNonZero++] = iColumn;
752	array[iColumn] = -value;
753	}
754	}
755	} else {
756	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
757	double value = `0.0`;
758	CoinBigIndex j;
759	const double * columnScale = model->columnScale();
760	for (j = columnStart[iColumn];
761	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
762	int iRow = row[j];
763	value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
764	}
765	value = scalar columnScale[iColumn];
766	if (fabs(value) > zeroTolerance) {
767	index[numberNonZero++] = iColumn;
768	array[iColumn] = value;
769	}
770	}
771	}
772	}
773	}
774	columnArray->setNumElements(numberNonZero);
775	y->setNumElements(`0`);
776	} else {
777	// do by row
778	rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
779	}
780	if (packed)
781	columnArray->setPackedMode(true);
782	if (`0`) {
783	columnArray->checkClean();
784	int numberNonZero = columnArray->getNumElements();
785	int * index = columnArray->getIndices();
786	double * array = columnArray->denseVector();
787	int i;
788	for (i = `0`; i < numberNonZero; i++) {
789	int j = index[i];
790	double value;
791	if (packed)
792	value = array[i];
793	else
794	value = array[j];
795	printf("Ti %d %d %g\n", i, j, value);
796	}
797	}
798	}
799	//static int xA=0;
800	//static int xB=0;
801	//static int xC=0;
802	//static int xD=0;
803	//static double yA=0.0;
804	//static double yC=0.0;
805	/ Return <code>x * scalar * A + y</code> in <code>z</code>.*
806	Note - If x packed mode - then z packed mode
807	This does by column and knows no gaps
808	Squashes small elements and knows about ClpSimplex /*
809	void
810	ClpPackedMatrix::transposeTimesByColumn(const ClpSimplex * model, double scalar,
811	const CoinIndexedVector * rowArray,
812	CoinIndexedVector * y,
813	CoinIndexedVector * columnArray) const
814	{
815	double * COIN_RESTRICT pi = rowArray->denseVector();
816	int numberNonZero = `0`;
817	int * COIN_RESTRICT index = columnArray->getIndices();
818	double * COIN_RESTRICT array = columnArray->denseVector();
819	int numberInRowArray = rowArray->getNumElements();
820	// maybe I need one in OsiSimplex
821	double zeroTolerance = model->zeroTolerance();
822	bool packed = rowArray->packedMode();
823	// do by column
824	int iColumn;
825	// get matrix data pointers
826	const int * COIN_RESTRICT row = matrix_->getIndices();
827	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
828	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
829	const double * COIN_RESTRICT rowScale = model->rowScale();
830	assert (!y->getNumElements());
831	assert (numberActiveColumns_ > `0`);
832	const ClpPackedMatrix * thisMatrix = this;
833	#if 0
834	ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
835	if (rowScale && scaledMatrix) {
836	rowScale = NULL;
837	// get matrix data pointers
838	row = scaledMatrix->getIndices();
839	columnStart = scaledMatrix->getVectorStarts();
840	elementByColumn = scaledMatrix->getElements();
841	thisMatrix = scaledMatrix;
842	//printf("scaledMatrix\n");
843	} else if (rowScale) {
844	//printf("no scaledMatrix\n");
845	} else {
846	//printf("no rowScale\n");
847	}
848	#endif
849	if (packed) {
850	// need to expand pi into y
851	assert(y->capacity() >= model->numberRows());
852	double * piOld = pi;
853	pi = y->denseVector();
854	const int * COIN_RESTRICT whichRow = rowArray->getIndices();
855	int i;
856	if (!rowScale) {
857	// modify pi so can collapse to one loop
858	if (scalar == -`1.0`) {
859	//yA += numberInRowArray;
860	for (i = `0`; i < numberInRowArray; i++) {
861	int iRow = whichRow[i];
862	pi[iRow] = -piOld[i];
863	}
864	} else {
865	for (i = `0`; i < numberInRowArray; i++) {
866	int iRow = whichRow[i];
867	pi[iRow] = scalar * piOld[i];
868	}
869	}
870	if (!columnCopy_) {
871	if ((model->specialOptions(), `131072`) != `0`) {
872	if(model->spareIntArray_[`0`] > `0`) {
873	CoinIndexedVector * spareArray = model->rowArray(`3`);
874	// also do dualColumn stuff
875	double * spare = spareArray->denseVector();
876	int * spareIndex = spareArray->getIndices();
877	const double * reducedCost = model->djRegion(`0`);
878	double multiplier[] = { -`1.0`, `1.0`};
879	double dualT = - model->currentDualTolerance();
880	double acceptablePivot = model->spareDoubleArray_[`0`];
881	// We can also see if infeasible or pivoting on free
882	double tentativeTheta = `1.0e15`;
883	double upperTheta = `1.0e31`;
884	double bestPossible = `0.0`;
885	int addSequence = model->numberColumns();
886	const unsigned char * statusArray = model->statusArray() + addSequence;
887	int numberRemaining = `0`;
888	assert (scalar == -`1.0`);
889	for (i = `0`; i < numberInRowArray; i++) {
890	int iSequence = whichRow[i];
891	int iStatus = (statusArray[iSequence] & `3`) - `1`;
892	if (iStatus) {
893	double mult = multiplier[iStatus-`1`];
894	double alpha = piOld[i] * mult;
895	double oldValue;
896	double value;
897	if (alpha > `0.0`) {
898	oldValue = reducedCost[iSequence] * mult;
899	value = oldValue - tentativeTheta * alpha;
900	if (value < dualT) {
901	bestPossible = CoinMax(bestPossible, alpha);
902	value = oldValue - upperTheta * alpha;
903	if (value < dualT && alpha >= acceptablePivot) {
904	upperTheta = (oldValue - dualT) / alpha;
905	//tentativeTheta = CoinMin(2.0upperTheta,tentativeTheta);*
906	}
907	// add to list
908	spare[numberRemaining] = alpha * mult;
909	spareIndex[numberRemaining++] = iSequence + addSequence;
910	}
911	}
912	}
913	}
914	numberNonZero =
915	thisMatrix->gutsOfTransposeTimesUnscaled(pi,
916	columnArray->getIndices(),
917	columnArray->denseVector(),
918	model->statusArray(),
919	spareIndex,
920	spare,
921	model->djRegion(`1`),
922	upperTheta,
923	bestPossible,
924	acceptablePivot,
925	model->currentDualTolerance(),
926	numberRemaining,
927	zeroTolerance);
928	model->spareDoubleArray_[`0`] = upperTheta;
929	model->spareDoubleArray_[`1`] = bestPossible;
930	spareArray->setNumElements(numberRemaining);
931	// signal partially done
932	model->spareIntArray_[`0`] = -`2`;
933	} else {
934	numberNonZero =
935	thisMatrix->gutsOfTransposeTimesUnscaled(pi,
936	columnArray->getIndices(),
937	columnArray->denseVector(),
938	model->statusArray(),
939	zeroTolerance);
940	}
941	} else {
942	numberNonZero =
943	thisMatrix->gutsOfTransposeTimesUnscaled(pi,
944	columnArray->getIndices(),
945	columnArray->denseVector(),
946	zeroTolerance);
947	}
948	columnArray->setNumElements(numberNonZero);
949	//xA++;
950	} else {
951	columnCopy_->transposeTimes(model, pi, columnArray);
952	numberNonZero = columnArray->getNumElements();
953	//xB++;
954	}
955	} else {
956	#ifdef CLP_INVESTIGATE
957	if (model->clpScaledMatrix())
958	printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
959	#endif
960	// scaled
961	// modify pi so can collapse to one loop
962	if (scalar == -`1.0`) {
963	//yC += numberInRowArray;
964	for (i = `0`; i < numberInRowArray; i++) {
965	int iRow = whichRow[i];
966	pi[iRow] = -piOld[i] * rowScale[iRow];
967	}
968	} else {
969	for (i = `0`; i < numberInRowArray; i++) {
970	int iRow = whichRow[i];
971	pi[iRow] = scalar * piOld[i] * rowScale[iRow];
972	}
973	}
974	const double * columnScale = model->columnScale();
975	if (!columnCopy_) {
976	if ((model->specialOptions(), `131072`) != `0`)
977	numberNonZero =
978	gutsOfTransposeTimesScaled(pi, columnScale,
979	columnArray->getIndices(),
980	columnArray->denseVector(),
981	model->statusArray(),
982	zeroTolerance);
983	else
984	numberNonZero =
985	gutsOfTransposeTimesScaled(pi, columnScale,
986	columnArray->getIndices(),
987	columnArray->denseVector(),
988	zeroTolerance);
989	columnArray->setNumElements(numberNonZero);
990	//xC++;
991	} else {
992	columnCopy_->transposeTimes(model, pi, columnArray);
993	numberNonZero = columnArray->getNumElements();
994	//xD++;
995	}
996	}
997	// zero out
998	int numberRows = model->numberRows();
999	if (numberInRowArray * `4` < numberRows) {
1000	for (i = `0`; i < numberInRowArray; i++) {
1001	int iRow = whichRow[i];
1002	pi[iRow] = `0.0`;
1003	}
1004	} else {
1005	CoinZeroN(pi, numberRows);
1006	}
1007	//int kA=xA+xB;
1008	//int kC=xC+xD;
1009	//if ((kA+kC)%10000==0)
1010	//printf("AA %d %d %g, CC %d %d %g\n",
1011	// xA,xB,kA ? yA/(double)(kA): 0.0,xC,xD,kC ? yC/(double) (kC) :0.0);
1012	} else {
1013	if (!rowScale) {
1014	if (scalar == -`1.0`) {
1015	double value = `0.0`;
1016	CoinBigIndex j;
1017	CoinBigIndex end = columnStart[`1`];
1018	for (j = columnStart[`0`]; j < end; j++) {
1019	int iRow = row[j];
1020	value += pi[iRow] * elementByColumn[j];
1021	}
1022	for (iColumn = `0`; iColumn < numberActiveColumns_ - `1`; iColumn++) {
1023	CoinBigIndex start = end;
1024	end = columnStart[iColumn+`2`];
1025	if (fabs(value) > zeroTolerance) {
1026	array[iColumn] = -value;
1027	index[numberNonZero++] = iColumn;
1028	}
1029	value = `0.0`;
1030	for (j = start; j < end; j++) {
1031	int iRow = row[j];
1032	value += pi[iRow] * elementByColumn[j];
1033	}
1034	}
1035	if (fabs(value) > zeroTolerance) {
1036	array[iColumn] = -value;
1037	index[numberNonZero++] = iColumn;
1038	}
1039	} else {
1040	double value = `0.0`;
1041	CoinBigIndex j;
1042	CoinBigIndex end = columnStart[`1`];
1043	for (j = columnStart[`0`]; j < end; j++) {
1044	int iRow = row[j];
1045	value += pi[iRow] * elementByColumn[j];
1046	}
1047	for (iColumn = `0`; iColumn < numberActiveColumns_ - `1`; iColumn++) {
1048	value *= scalar;
1049	CoinBigIndex start = end;
1050	end = columnStart[iColumn+`2`];
1051	if (fabs(value) > zeroTolerance) {
1052	array[iColumn] = value;
1053	index[numberNonZero++] = iColumn;
1054	}
1055	value = `0.0`;
1056	for (j = start; j < end; j++) {
1057	int iRow = row[j];
1058	value += pi[iRow] * elementByColumn[j];
1059	}
1060	}
1061	value *= scalar;
1062	if (fabs(value) > zeroTolerance) {
1063	array[iColumn] = value;
1064	index[numberNonZero++] = iColumn;
1065	}
1066	}
1067	} else {
1068	#ifdef CLP_INVESTIGATE
1069	if (model->clpScaledMatrix())
1070	printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
1071	#endif
1072	// scaled
1073	if (scalar == -`1.0`) {
1074	const double * columnScale = model->columnScale();
1075	double value = `0.0`;
1076	double scale = columnScale[`0`];
1077	CoinBigIndex j;
1078	CoinBigIndex end = columnStart[`1`];
1079	for (j = columnStart[`0`]; j < end; j++) {
1080	int iRow = row[j];
1081	value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1082	}
1083	for (iColumn = `0`; iColumn < numberActiveColumns_ - `1`; iColumn++) {
1084	value *= scale;
1085	CoinBigIndex start = end;
1086	end = columnStart[iColumn+`2`];
1087	scale = columnScale[iColumn+`1`];
1088	if (fabs(value) > zeroTolerance) {
1089	array[iColumn] = -value;
1090	index[numberNonZero++] = iColumn;
1091	}
1092	value = `0.0`;
1093	for (j = start; j < end; j++) {
1094	int iRow = row[j];
1095	value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1096	}
1097	}
1098	value *= scale;
1099	if (fabs(value) > zeroTolerance) {
1100	array[iColumn] = -value;
1101	index[numberNonZero++] = iColumn;
1102	}
1103	} else {
1104	const double * columnScale = model->columnScale();
1105	double value = `0.0`;
1106	double scale = columnScale[`0`] * scalar;
1107	CoinBigIndex j;
1108	CoinBigIndex end = columnStart[`1`];
1109	for (j = columnStart[`0`]; j < end; j++) {
1110	int iRow = row[j];
1111	value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1112	}
1113	for (iColumn = `0`; iColumn < numberActiveColumns_ - `1`; iColumn++) {
1114	value *= scale;
1115	CoinBigIndex start = end;
1116	end = columnStart[iColumn+`2`];
1117	scale = columnScale[iColumn+`1`] * scalar;
1118	if (fabs(value) > zeroTolerance) {
1119	array[iColumn] = value;
1120	index[numberNonZero++] = iColumn;
1121	}
1122	value = `0.0`;
1123	for (j = start; j < end; j++) {
1124	int iRow = row[j];
1125	value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1126	}
1127	}
1128	value *= scale;
1129	if (fabs(value) > zeroTolerance) {
1130	array[iColumn] = value;
1131	index[numberNonZero++] = iColumn;
1132	}
1133	}
1134	}
1135	}
1136	columnArray->setNumElements(numberNonZero);
1137	y->setNumElements(`0`);
1138	if (packed)
1139	columnArray->setPackedMode(true);
1140	}
1141	/ Return <code>x * A + y</code> in <code>z</code>.*
1142	Squashes small elements and knows about ClpSimplex /*
1143	void
1144	ClpPackedMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
1145	const CoinIndexedVector * rowArray,
1146	CoinIndexedVector * y,
1147	CoinIndexedVector * columnArray) const
1148	{
1149	columnArray->clear();
1150	double * pi = rowArray->denseVector();
1151	int numberNonZero = `0`;
1152	int * index = columnArray->getIndices();
1153	double * array = columnArray->denseVector();
1154	int numberInRowArray = rowArray->getNumElements();
1155	// maybe I need one in OsiSimplex
1156	double zeroTolerance = model->zeroTolerance();
1157	const int * column = matrix_->getIndices();
1158	const CoinBigIndex * rowStart = getVectorStarts();
1159	const double * element = getElements();
1160	const int * whichRow = rowArray->getIndices();
1161	bool packed = rowArray->packedMode();
1162	if (numberInRowArray > `2`) {
1163	// do by rows
1164	// * Row copy is already scaled*
1165	int iRow;
1166	int i;
1167	int numberOriginal = `0`;
1168	if (packed) {
1169	int * index = columnArray->getIndices();
1170	double * array = columnArray->denseVector();
1171	#if 0
1172	{
1173	double * array2 = y->denseVector();
1174	int numberColumns = matrix_->getNumCols();
1175	for (int i=`0`;i<numberColumns;i++) {
1176	assert(!array[i]);
1177	assert(!array2[i]);
1178	}
1179	}
1180	#endif
1181	//#define COIN_SPARSE_MATRIX 1
1182	#if COIN_SPARSE_MATRIX
1183	assert (!y->getNumElements());
1184	#if COIN_SPARSE_MATRIX != 2
1185	// and set up mark as char array
1186	char * marked = reinterpret_cast<char *> (index+columnArray->capacity());
1187	int * lookup = y->getIndices();
1188	#ifndef NDEBUG
1189	//int numberColumns = matrix_->getNumCols();
1190	//for (int i=0;i<numberColumns;i++)
1191	//assert(!marked[i]);
1192	#endif
1193	numberNonZero=gutsOfTransposeTimesByRowGE3a(rowArray,index,array,
1194	lookup,marked,zeroTolerance,scalar);
1195	#else
1196	double * array2 = y->denseVector();
1197	numberNonZero=gutsOfTransposeTimesByRowGE3(rowArray,index,array,
1198	array2,zeroTolerance,scalar);
1199	#endif
1200	#else
1201	int numberColumns = matrix_->getNumCols();
1202	numberNonZero = gutsOfTransposeTimesByRowGEK(rowArray, index, array,
1203	numberColumns, zeroTolerance, scalar);
1204	#endif
1205	columnArray->setNumElements(numberNonZero);
1206	} else {
1207	double * markVector = y->denseVector();
1208	numberNonZero = `0`;
1209	// and set up mark as char array
1210	char * marked = reinterpret_cast<char *> (markVector);
1211	for (i = `0`; i < numberOriginal; i++) {
1212	int iColumn = index[i];
1213	marked[iColumn] = `0`;
1214	}
1215
1216	for (i = `0`; i < numberInRowArray; i++) {
1217	iRow = whichRow[i];
1218	double value = pi[iRow] * scalar;
1219	CoinBigIndex j;
1220	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
1221	int iColumn = column[j];
1222	if (!marked[iColumn]) {
1223	marked[iColumn] = `1`;
1224	index[numberNonZero++] = iColumn;
1225	}
1226	array[iColumn] += value * element[j];
1227	}
1228	}
1229	// get rid of tiny values and zero out marked
1230	numberOriginal = numberNonZero;
1231	numberNonZero = `0`;
1232	for (i = `0`; i < numberOriginal; i++) {
1233	int iColumn = index[i];
1234	marked[iColumn] = `0`;
1235	if (fabs(array[iColumn]) > zeroTolerance) {
1236	index[numberNonZero++] = iColumn;
1237	} else {
1238	array[iColumn] = `0.0`;
1239	}
1240	}
1241	}
1242	} else if (numberInRowArray == `2`) {
1243	// do by rows when two rows
1244	int numberOriginal;
1245	int i;
1246	CoinBigIndex j;
1247	numberNonZero = `0`;
1248
1249	double value;
1250	if (packed) {
1251	gutsOfTransposeTimesByRowEQ2(rowArray, columnArray, y, zeroTolerance, scalar);
1252	numberNonZero = columnArray->getNumElements();
1253	} else {
1254	int iRow = whichRow[`0`];
1255	value = pi[iRow] * scalar;
1256	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
1257	int iColumn = column[j];
1258	double value2 = value * element[j];
1259	index[numberNonZero++] = iColumn;
1260	array[iColumn] = value2;
1261	}
1262	iRow = whichRow[`1`];
1263	value = pi[iRow] * scalar;
1264	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
1265	int iColumn = column[j];
1266	double value2 = value * element[j];
1267	// I am assuming no zeros in matrix
1268	if (array[iColumn])
1269	value2 += array[iColumn];
1270	else
1271	index[numberNonZero++] = iColumn;
1272	array[iColumn] = value2;
1273	}
1274	// get rid of tiny values and zero out marked
1275	numberOriginal = numberNonZero;
1276	numberNonZero = `0`;
1277	for (i = `0`; i < numberOriginal; i++) {
1278	int iColumn = index[i];
1279	if (fabs(array[iColumn]) > zeroTolerance) {
1280	index[numberNonZero++] = iColumn;
1281	} else {
1282	array[iColumn] = `0.0`;
1283	}
1284	}
1285	}
1286	} else if (numberInRowArray == `1`) {
1287	// Just one row
1288	int iRow = rowArray->getIndices()[`0`];
1289	numberNonZero = `0`;
1290	CoinBigIndex j;
1291	if (packed) {
1292	gutsOfTransposeTimesByRowEQ1(rowArray, columnArray, zeroTolerance, scalar);
1293	numberNonZero = columnArray->getNumElements();
1294	} else {
1295	double value = pi[iRow] * scalar;
1296	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
1297	int iColumn = column[j];
1298	double value2 = value * element[j];
1299	if (fabs(value2) > zeroTolerance) {
1300	index[numberNonZero++] = iColumn;
1301	array[iColumn] = value2;
1302	}
1303	}
1304	}
1305	}
1306	columnArray->setNumElements(numberNonZero);
1307	y->setNumElements(`0`);
1308	}
1309	// Meat of transposeTimes by column when not scaled
1310	int
1311	ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
1312	int * COIN_RESTRICT index,
1313	double * COIN_RESTRICT array,
1314	const double zeroTolerance) const
1315	{
1316	int numberNonZero = `0`;
1317	// get matrix data pointers
1318	const int * COIN_RESTRICT row = matrix_->getIndices();
1319	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1320	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1321	#if 1 //ndef INTEL_MKL
1322	double value = `0.0`;
1323	CoinBigIndex j;
1324	CoinBigIndex end = columnStart[`1`];
1325	for (j = columnStart[`0`]; j < end; j++) {
1326	int iRow = row[j];
1327	value += pi[iRow] * elementByColumn[j];
1328	}
1329	int iColumn;
1330	for (iColumn = `0`; iColumn < numberActiveColumns_ - `1`; iColumn++) {
1331	CoinBigIndex start = end;
1332	end = columnStart[iColumn+`2`];
1333	if (fabs(value) > zeroTolerance) {
1334	array[numberNonZero] = value;
1335	index[numberNonZero++] = iColumn;
1336	}
1337	value = `0.0`;
1338	for (j = start; j < end; j++) {
1339	int iRow = row[j];
1340	value += pi[iRow] * elementByColumn[j];
1341	}
1342	}
1343	if (fabs(value) > zeroTolerance) {
1344	array[numberNonZero] = value;
1345	index[numberNonZero++] = iColumn;
1346	}
1347	#else
1348	char transA = `'N'`;
1349	//int numberRows = matrix_->getNumRows();
1350	mkl_cspblas_dcsrgemv(&transA, const_cast<int *>(&numberActiveColumns_),
1351	const_cast<double *>(elementByColumn),
1352	const_cast<int *>(columnStart),
1353	const_cast<int *>(row),
1354	const_cast<double *>(pi), array);
1355	int iColumn;
1356	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
1357	double value = array[iColumn];
1358	if (value) {
1359	array[iColumn] = `0.0`;
1360	if (fabs(value) > zeroTolerance) {
1361	array[numberNonZero] = value;
1362	index[numberNonZero++] = iColumn;
1363	}
1364	}
1365	}
1366	#endif
1367	return numberNonZero;
1368	}
1369	// Meat of transposeTimes by column when scaled
1370	int
1371	ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,
1372	const double * COIN_RESTRICT columnScale,
1373	int * COIN_RESTRICT index,
1374	double * COIN_RESTRICT array,
1375	const double zeroTolerance) const
1376	{
1377	int numberNonZero = `0`;
1378	// get matrix data pointers
1379	const int * COIN_RESTRICT row = matrix_->getIndices();
1380	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1381	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1382	double value = `0.0`;
1383	double scale = columnScale[`0`];
1384	CoinBigIndex j;
1385	CoinBigIndex end = columnStart[`1`];
1386	for (j = columnStart[`0`]; j < end; j++) {
1387	int iRow = row[j];
1388	value += pi[iRow] * elementByColumn[j];
1389	}
1390	int iColumn;
1391	for (iColumn = `0`; iColumn < numberActiveColumns_ - `1`; iColumn++) {
1392	value *= scale;
1393	CoinBigIndex start = end;
1394	scale = columnScale[iColumn+`1`];
1395	end = columnStart[iColumn+`2`];
1396	if (fabs(value) > zeroTolerance) {
1397	array[numberNonZero] = value;
1398	index[numberNonZero++] = iColumn;
1399	}
1400	value = `0.0`;
1401	for (j = start; j < end; j++) {
1402	int iRow = row[j];
1403	value += pi[iRow] * elementByColumn[j];
1404	}
1405	}
1406	value *= scale;
1407	if (fabs(value) > zeroTolerance) {
1408	array[numberNonZero] = value;
1409	index[numberNonZero++] = iColumn;
1410	}
1411	return numberNonZero;
1412	}
1413	// Meat of transposeTimes by column when not scaled
1414	int
1415	ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
1416	int * COIN_RESTRICT index,
1417	double * COIN_RESTRICT array,
1418	const unsigned char * COIN_RESTRICT status,
1419	const double zeroTolerance) const
1420	{
1421	int numberNonZero = `0`;
1422	// get matrix data pointers
1423	const int * COIN_RESTRICT row = matrix_->getIndices();
1424	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1425	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1426	double value = `0.0`;
1427	int jColumn = -`1`;
1428	for (int iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
1429	bool wanted = ((status[iColumn] & `3`) != `1`);
1430	if (fabs(value) > zeroTolerance) {
1431	array[numberNonZero] = value;
1432	index[numberNonZero++] = jColumn;
1433	}
1434	value = `0.0`;
1435	if (wanted) {
1436	CoinBigIndex start = columnStart[iColumn];
1437	CoinBigIndex end = columnStart[iColumn+`1`];
1438	jColumn = iColumn;
1439	int n = end - start;
1440	bool odd = (n & `1`) != `0`;
1441	n = n >> `1`;
1442	const int * COIN_RESTRICT rowThis = row + start;
1443	const double * COIN_RESTRICT elementThis = elementByColumn + start;
1444	for (; n; n--) {
1445	int iRow0 = *rowThis;
1446	int iRow1 = *(rowThis + `1`);
1447	rowThis += `2`;
1448	value += pi[iRow0] * (*elementThis);
1449	value += pi[iRow1] * (*(elementThis + `1`));
1450	elementThis += `2`;
1451	}
1452	if (odd) {
1453	int iRow = *rowThis;
1454	value += pi[iRow] * (*elementThis);
1455	}
1456	}
1457	}
1458	if (fabs(value) > zeroTolerance) {
1459	array[numberNonZero] = value;
1460	index[numberNonZero++] = jColumn;
1461	}
1462	return numberNonZero;
1463	}
1464	/ Meat of transposeTimes by column when not scaled and skipping*
1465	and doing part of dualColumn /*
1466	int
1467	ClpPackedMatrix::gutsOfTransposeTimesUnscaled(const double * COIN_RESTRICT pi,
1468	int * COIN_RESTRICT index,
1469	double * COIN_RESTRICT array,
1470	const unsigned char * COIN_RESTRICT status,
1471	int * COIN_RESTRICT spareIndex,
1472	double * COIN_RESTRICT spareArray,
1473	const double * COIN_RESTRICT reducedCost,
1474	double & upperThetaP,
1475	double & bestPossibleP,
1476	double acceptablePivot,
1477	double dualTolerance,
1478	int & numberRemainingP,
1479	const double zeroTolerance) const
1480	{
1481	double tentativeTheta = `1.0e15`;
1482	int numberRemaining = numberRemainingP;
1483	double upperTheta = upperThetaP;
1484	double bestPossible = bestPossibleP;
1485	int numberNonZero = `0`;
1486	// get matrix data pointers
1487	const int * COIN_RESTRICT row = matrix_->getIndices();
1488	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1489	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1490	double multiplier[] = { -`1.0`, `1.0`};
1491	double dualT = - dualTolerance;
1492	for (int iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
1493	int wanted = (status[iColumn] & `3`) - `1`;
1494	if (wanted) {
1495	double value = `0.0`;
1496	CoinBigIndex start = columnStart[iColumn];
1497	CoinBigIndex end = columnStart[iColumn+`1`];
1498	int n = end - start;
1499	#if 1
1500	bool odd = (n & `1`) != `0`;
1501	n = n >> `1`;
1502	const int * COIN_RESTRICT rowThis = row + start;
1503	const double * COIN_RESTRICT elementThis = elementByColumn + start;
1504	for (; n; n--) {
1505	int iRow0 = *rowThis;
1506	int iRow1 = *(rowThis + `1`);
1507	rowThis += `2`;
1508	value += pi[iRow0] * (*elementThis);
1509	value += pi[iRow1] * (*(elementThis + `1`));
1510	elementThis += `2`;
1511	}
1512	if (odd) {
1513	int iRow = *rowThis;
1514	value += pi[iRow] * (*elementThis);
1515	}
1516	#else
1517	const int * COIN_RESTRICT rowThis = &row[end-`16`];
1518	const double * COIN_RESTRICT elementThis = &elementByColumn[end-`16`];
1519	bool odd = (n & `1`) != `0`;
1520	n = n >> `1`;
1521	double value2 = `0.0`;
1522	if (odd) {
1523	int iRow = row[start];
1524	value2 = pi[iRow] * elementByColumn[start];
1525	}
1526	switch (n) {
1527	default: {
1528	if (odd) {
1529	start++;
1530	}
1531	n -= `8`;
1532	for (; n; n--) {
1533	int iRow0 = row[start];
1534	int iRow1 = row[start+`1`];
1535	value += pi[iRow0] * elementByColumn[start];
1536	value2 += pi[iRow1] * elementByColumn[start+`1`];
1537	start += `2`;
1538	}
1539	case `8`: {
1540	int iRow0 = rowThis[`16`-`16`];
1541	int iRow1 = rowThis[`16`-`15`];
1542	value += pi[iRow0] * elementThis[`16`-`16`];
1543	value2 += pi[iRow1] * elementThis[`16`-`15`];
1544	}
1545	case `7`: {
1546	int iRow0 = rowThis[`16`-`14`];
1547	int iRow1 = rowThis[`16`-`13`];
1548	value += pi[iRow0] * elementThis[`16`-`14`];
1549	value2 += pi[iRow1] * elementThis[`16`-`13`];
1550	}
1551	case `6`: {
1552	int iRow0 = rowThis[`16`-`12`];
1553	int iRow1 = rowThis[`16`-`11`];
1554	value += pi[iRow0] * elementThis[`16`-`12`];
1555	value2 += pi[iRow1] * elementThis[`16`-`11`];
1556	}
1557	case `5`: {
1558	int iRow0 = rowThis[`16`-`10`];
1559	int iRow1 = rowThis[`16`-`9`];
1560	value += pi[iRow0] * elementThis[`16`-`10`];
1561	value2 += pi[iRow1] * elementThis[`16`-`9`];
1562	}
1563	case `4`: {
1564	int iRow0 = rowThis[`16`-`8`];
1565	int iRow1 = rowThis[`16`-`7`];
1566	value += pi[iRow0] * elementThis[`16`-`8`];
1567	value2 += pi[iRow1] * elementThis[`16`-`7`];
1568	}
1569	case `3`: {
1570	int iRow0 = rowThis[`16`-`6`];
1571	int iRow1 = rowThis[`16`-`5`];
1572	value += pi[iRow0] * elementThis[`16`-`6`];
1573	value2 += pi[iRow1] * elementThis[`16`-`5`];
1574	}
1575	case `2`: {
1576	int iRow0 = rowThis[`16`-`4`];
1577	int iRow1 = rowThis[`16`-`3`];
1578	value += pi[iRow0] * elementThis[`16`-`4`];
1579	value2 += pi[iRow1] * elementThis[`16`-`3`];
1580	}
1581	case `1`: {
1582	int iRow0 = rowThis[`16`-`2`];
1583	int iRow1 = rowThis[`16`-`1`];
1584	value += pi[iRow0] * elementThis[`16`-`2`];
1585	value2 += pi[iRow1] * elementThis[`16`-`1`];
1586	}
1587	case `0`:
1588	;
1589	}
1590	}
1591	value += value2;
1592	#endif
1593	if (fabs(value) > zeroTolerance) {
1594	double mult = multiplier[wanted-`1`];
1595	double alpha = value * mult;
1596	array[numberNonZero] = value;
1597	index[numberNonZero++] = iColumn;
1598	if (alpha > `0.0`) {
1599	double oldValue = reducedCost[iColumn] * mult;
1600	double value = oldValue - tentativeTheta * alpha;
1601	if (value < dualT) {
1602	bestPossible = CoinMax(bestPossible, alpha);
1603	value = oldValue - upperTheta * alpha;
1604	if (value < dualT && alpha >= acceptablePivot) {
1605	upperTheta = (oldValue - dualT) / alpha;
1606	//tentativeTheta = CoinMin(2.0upperTheta,tentativeTheta);*
1607	}
1608	// add to list
1609	spareArray[numberRemaining] = alpha * mult;
1610	spareIndex[numberRemaining++] = iColumn;
1611	}
1612	}
1613	}
1614	}
1615	}
1616	numberRemainingP = numberRemaining;
1617	upperThetaP = upperTheta;
1618	bestPossibleP = bestPossible;
1619	return numberNonZero;
1620	}
1621	// Meat of transposeTimes by column when scaled
1622	int
1623	ClpPackedMatrix::gutsOfTransposeTimesScaled(const double * COIN_RESTRICT pi,
1624	const double * COIN_RESTRICT columnScale,
1625	int * COIN_RESTRICT index,
1626	double * COIN_RESTRICT array,
1627	const unsigned char * COIN_RESTRICT status, const double zeroTolerance) const
1628	{
1629	int numberNonZero = `0`;
1630	// get matrix data pointers
1631	const int * COIN_RESTRICT row = matrix_->getIndices();
1632	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1633	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1634	double value = `0.0`;
1635	int jColumn = -`1`;
1636	for (int iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
1637	bool wanted = ((status[iColumn] & `3`) != `1`);
1638	if (fabs(value) > zeroTolerance) {
1639	array[numberNonZero] = value;
1640	index[numberNonZero++] = jColumn;
1641	}
1642	value = `0.0`;
1643	if (wanted) {
1644	double scale = columnScale[iColumn];
1645	CoinBigIndex start = columnStart[iColumn];
1646	CoinBigIndex end = columnStart[iColumn+`1`];
1647	jColumn = iColumn;
1648	for (CoinBigIndex j = start; j < end; j++) {
1649	int iRow = row[j];
1650	value += pi[iRow] * elementByColumn[j];
1651	}
1652	value *= scale;
1653	}
1654	}
1655	if (fabs(value) > zeroTolerance) {
1656	array[numberNonZero] = value;
1657	index[numberNonZero++] = jColumn;
1658	}
1659	return numberNonZero;
1660	}
1661	// Meat of transposeTimes by row n > K if packed - returns number nonzero
1662	int
1663	ClpPackedMatrix::gutsOfTransposeTimesByRowGEK(const CoinIndexedVector * COIN_RESTRICT piVector,
1664	int * COIN_RESTRICT index,
1665	double * COIN_RESTRICT output,
1666	int numberColumns,
1667	const double tolerance,
1668	const double scalar) const
1669	{
1670	const double * COIN_RESTRICT pi = piVector->denseVector();
1671	int numberInRowArray = piVector->getNumElements();
1672	const int * COIN_RESTRICT column = matrix_->getIndices();
1673	const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
1674	const double * COIN_RESTRICT element = matrix_->getElements();
1675	const int * COIN_RESTRICT whichRow = piVector->getIndices();
1676	// * Row copy is already scaled*
1677	for (int i = `0`; i < numberInRowArray; i++) {
1678	int iRow = whichRow[i];
1679	double value = pi[i] * scalar;
1680	CoinBigIndex start = rowStart[iRow];
1681	CoinBigIndex end = rowStart[iRow+`1`];
1682	int n = end - start;
1683	const int * COIN_RESTRICT columnThis = column + start;
1684	const double * COIN_RESTRICT elementThis = element + start;
1685
1686	// could do by twos
1687	for (; n; n--) {
1688	int iColumn = *columnThis;
1689	columnThis++;
1690	double elValue = *elementThis;
1691	elementThis++;
1692	elValue *= value;
1693	output[iColumn] += elValue;
1694	}
1695	}
1696	// get rid of tiny values and count
1697	int numberNonZero = `0`;
1698	for (int i = `0`; i < numberColumns; i++) {
1699	double value = output[i];
1700	if (value) {
1701	output[i] = `0.0`;
1702	if (fabs(value) > tolerance) {
1703	output[numberNonZero] = value;
1704	index[numberNonZero++] = i;
1705	}
1706	}
1707	}
1708	#ifndef NDEBUG
1709	for (int i = numberNonZero; i < numberColumns; i++)
1710	assert(!output[i]);
1711	#endif
1712	return numberNonZero;
1713	}
1714	// Meat of transposeTimes by row n == 2 if packed
1715	void
1716	ClpPackedMatrix::gutsOfTransposeTimesByRowEQ2(const CoinIndexedVector * piVector, CoinIndexedVector * output,
1717	CoinIndexedVector * spareVector, const double tolerance, const double scalar) const
1718	{
1719	double * pi = piVector->denseVector();
1720	int numberNonZero = `0`;
1721	int * index = output->getIndices();
1722	double * array = output->denseVector();
1723	const int * column = matrix_->getIndices();
1724	const CoinBigIndex * rowStart = matrix_->getVectorStarts();
1725	const double * element = matrix_->getElements();
1726	const int * whichRow = piVector->getIndices();
1727	int iRow0 = whichRow[`0`];
1728	int iRow1 = whichRow[`1`];
1729	double pi0 = pi[`0`];
1730	double pi1 = pi[`1`];
1731	if (rowStart[iRow0+`1`] - rowStart[iRow0] >
1732	rowStart[iRow1+`1`] - rowStart[iRow1]) {
1733	// do one with fewer first
1734	iRow0 = iRow1;
1735	iRow1 = whichRow[`0`];
1736	pi0 = pi1;
1737	pi1 = pi[`0`];
1738	}
1739	// and set up mark as char array
1740	char * marked = reinterpret_cast<char *> (index + output->capacity());
1741	int * lookup = spareVector->getIndices();
1742	double value = pi0 * scalar;
1743	CoinBigIndex j;
1744	for (j = rowStart[iRow0]; j < rowStart[iRow0+`1`]; j++) {
1745	int iColumn = column[j];
1746	double elValue = element[j];
1747	double value2 = value * elValue;
1748	array[numberNonZero] = value2;
1749	marked[iColumn] = `1`;
1750	lookup[iColumn] = numberNonZero;
1751	index[numberNonZero++] = iColumn;
1752	}
1753	int numberOriginal = numberNonZero;
1754	value = pi1 * scalar;
1755	for (j = rowStart[iRow1]; j < rowStart[iRow1+`1`]; j++) {
1756	int iColumn = column[j];
1757	double elValue = element[j];
1758	double value2 = value * elValue;
1759	// I am assuming no zeros in matrix
1760	if (marked[iColumn]) {
1761	int iLookup = lookup[iColumn];
1762	array[iLookup] += value2;
1763	} else {
1764	if (fabs(value2) > tolerance) {
1765	array[numberNonZero] = value2;
1766	index[numberNonZero++] = iColumn;
1767	}
1768	}
1769	}
1770	// get rid of tiny values and zero out marked
1771	int i;
1772	int iFirst = numberNonZero;
1773	for (i = `0`; i < numberOriginal; i++) {
1774	int iColumn = index[i];
1775	marked[iColumn] = `0`;
1776	if (fabs(array[i]) <= tolerance) {
1777	if (numberNonZero > numberOriginal) {
1778	numberNonZero--;
1779	double value = array[numberNonZero];
1780	array[numberNonZero] = `0.0`;
1781	array[i] = value;
1782	index[i] = index[numberNonZero];
1783	} else {
1784	iFirst = i;
1785	}
1786	}
1787	}
1788
1789	if (iFirst < numberNonZero) {
1790	int n = iFirst;
1791	for (i = n; i < numberOriginal; i++) {
1792	int iColumn = index[i];
1793	double value = array[i];
1794	array[i] = `0.0`;
1795	if (fabs(value) > tolerance) {
1796	array[n] = value;
1797	index[n++] = iColumn;
1798	}
1799	}
1800	for (; i < numberNonZero; i++) {
1801	int iColumn = index[i];
1802	double value = array[i];
1803	array[i] = `0.0`;
1804	array[n] = value;
1805	index[n++] = iColumn;
1806	}
1807	numberNonZero = n;
1808	}
1809	output->setNumElements(numberNonZero);
1810	spareVector->setNumElements(`0`);
1811	}
1812	// Meat of transposeTimes by row n == 1 if packed
1813	void
1814	ClpPackedMatrix::gutsOfTransposeTimesByRowEQ1(const CoinIndexedVector * piVector, CoinIndexedVector * output,
1815	const double tolerance, const double scalar) const
1816	{
1817	double * pi = piVector->denseVector();
1818	int numberNonZero = `0`;
1819	int * index = output->getIndices();
1820	double * array = output->denseVector();
1821	const int * column = matrix_->getIndices();
1822	const CoinBigIndex * rowStart = matrix_->getVectorStarts();
1823	const double * element = matrix_->getElements();
1824	int iRow = piVector->getIndices()[`0`];
1825	numberNonZero = `0`;
1826	CoinBigIndex j;
1827	double value = pi[`0`] * scalar;
1828	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
1829	int iColumn = column[j];
1830	double elValue = element[j];
1831	double value2 = value * elValue;
1832	if (fabs(value2) > tolerance) {
1833	array[numberNonZero] = value2;
1834	index[numberNonZero++] = iColumn;
1835	}
1836	}
1837	output->setNumElements(numberNonZero);
1838	}
1839	/ Return <code>x A in <code>z</code> but
1840	just for indices in y.
1841	Squashes small elements and knows about ClpSimplex /*
1842	void
1843	ClpPackedMatrix::subsetTransposeTimes(const ClpSimplex * model,
1844	const CoinIndexedVector * rowArray,
1845	const CoinIndexedVector * y,
1846	CoinIndexedVector * columnArray) const
1847	{
1848	columnArray->clear();
1849	double * COIN_RESTRICT pi = rowArray->denseVector();
1850	double * COIN_RESTRICT array = columnArray->denseVector();
1851	int jColumn;
1852	// get matrix data pointers
1853	const int * COIN_RESTRICT row = matrix_->getIndices();
1854	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
1855	const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
1856	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
1857	const double * COIN_RESTRICT rowScale = model->rowScale();
1858	int numberToDo = y->getNumElements();
1859	const int * COIN_RESTRICT which = y->getIndices();
1860	assert (!rowArray->packedMode());
1861	columnArray->setPacked();
1862	ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
1863	int flags = flags_;
1864	if (rowScale && scaledMatrix && !(scaledMatrix->flags() & `2`)) {
1865	flags = `0`;
1866	rowScale = NULL;
1867	// get matrix data pointers
1868	row = scaledMatrix->getIndices();
1869	columnStart = scaledMatrix->getVectorStarts();
1870	elementByColumn = scaledMatrix->getElements();
1871	}
1872	if (!(flags & `2`) && numberToDo>`2`) {
1873	// no gaps
1874	if (!rowScale) {
1875	int iColumn = which[`0`];
1876	double value = `0.0`;
1877	CoinBigIndex j;
1878	int columnNext = which[`1`];
1879	CoinBigIndex startNext=columnStart[columnNext];
1880	//coin_prefetch_const(row+startNext);
1881	//coin_prefetch_const(elementByColumn+startNext);
1882	CoinBigIndex endNext=columnStart[columnNext+`1`];
1883	for (j = columnStart[iColumn];
1884	j < columnStart[iColumn+`1`]; j++) {
1885	int iRow = row[j];
1886	value += pi[iRow] * elementByColumn[j];
1887	}
1888	for (jColumn = `0`; jColumn < numberToDo - `2`; jColumn++) {
1889	CoinBigIndex start = startNext;
1890	CoinBigIndex end = endNext;
1891	columnNext = which[jColumn+`2`];
1892	startNext=columnStart[columnNext];
1893	//coin_prefetch_const(row+startNext);
1894	//coin_prefetch_const(elementByColumn+startNext);
1895	endNext=columnStart[columnNext+`1`];
1896	array[jColumn] = value;
1897	value = `0.0`;
1898	for (j = start; j < end; j++) {
1899	int iRow = row[j];
1900	value += pi[iRow] * elementByColumn[j];
1901	}
1902	}
1903	array[jColumn++] = value;
1904	value = `0.0`;
1905	for (j = startNext; j < endNext; j++) {
1906	int iRow = row[j];
1907	value += pi[iRow] * elementByColumn[j];
1908	}
1909	array[jColumn] = value;
1910	} else {
1911	#ifdef CLP_INVESTIGATE
1912	if (model->clpScaledMatrix())
1913	printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
1914	#endif
1915	// scaled
1916	const double * columnScale = model->columnScale();
1917	int iColumn = which[`0`];
1918	double value = `0.0`;
1919	double scale = columnScale[iColumn];
1920	CoinBigIndex j;
1921	for (j = columnStart[iColumn];
1922	j < columnStart[iColumn+`1`]; j++) {
1923	int iRow = row[j];
1924	value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1925	}
1926	for (jColumn = `0`; jColumn < numberToDo - `1`; jColumn++) {
1927	int iColumn = which[jColumn+`1`];
1928	value *= scale;
1929	scale = columnScale[iColumn];
1930	CoinBigIndex start = columnStart[iColumn];
1931	CoinBigIndex end = columnStart[iColumn+`1`];
1932	array[jColumn] = value;
1933	value = `0.0`;
1934	for (j = start; j < end; j++) {
1935	int iRow = row[j];
1936	value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1937	}
1938	}
1939	value *= scale;
1940	array[jColumn] = value;
1941	}
1942	} else if (numberToDo) {
1943	// gaps
1944	if (!rowScale) {
1945	for (jColumn = `0`; jColumn < numberToDo; jColumn++) {
1946	int iColumn = which[jColumn];
1947	double value = `0.0`;
1948	CoinBigIndex j;
1949	for (j = columnStart[iColumn];
1950	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1951	int iRow = row[j];
1952	value += pi[iRow] * elementByColumn[j];
1953	}
1954	array[jColumn] = value;
1955	}
1956	} else {
1957	#ifdef CLP_INVESTIGATE
1958	if (model->clpScaledMatrix())
1959	printf("scaledMatrix_ at %d of ClpPackedMatrix - flags %d (%d) n %d\n",
1960	__LINE__, flags_, model->clpScaledMatrix()->flags(), numberToDo);
1961	#endif
1962	// scaled
1963	const double * columnScale = model->columnScale();
1964	for (jColumn = `0`; jColumn < numberToDo; jColumn++) {
1965	int iColumn = which[jColumn];
1966	double value = `0.0`;
1967	CoinBigIndex j;
1968	for (j = columnStart[iColumn];
1969	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
1970	int iRow = row[j];
1971	value += pi[iRow] * elementByColumn[j] * rowScale[iRow];
1972	}
1973	value *= columnScale[iColumn];
1974	array[jColumn] = value;
1975	}
1976	}
1977	}
1978	}
1979	/ Returns true if can combine transposeTimes and subsetTransposeTimes*
1980	and if it would be faster /*
1981	bool
1982	ClpPackedMatrix::canCombine(const ClpSimplex * model,
1983	const CoinIndexedVector * pi) const
1984	{
1985	int numberInRowArray = pi->getNumElements();
1986	int numberRows = model->numberRows();
1987	bool packed = pi->packedMode();
1988	// factor should be smaller if doing both with two pi vectors
1989	double factor = `0.30`;
1990	// We may not want to do by row if there may be cache problems
1991	// It would be nice to find L2 cache size - for moment 512K
1992	// Be slightly optimistic
1993	if (numberActiveColumns_ * sizeof(double) > `1000000`) {
1994	if (numberRows * `10` < numberActiveColumns_)
1995	factor *= `0.333333333`;
1996	else if (numberRows * `4` < numberActiveColumns_)
1997	factor *= `0.5`;
1998	else if (numberRows * `2` < numberActiveColumns_)
1999	factor *= `0.66666666667`;
2000	//if (model->numberIterations()%50==0)
2001	//printf("%d nonzero\n",numberInRowArray);
2002	}
2003	// if not packed then bias a bit more towards by column
2004	if (!packed)
2005	factor *= `0.9`;
2006	return ((numberInRowArray > factor * numberRows \|\| !model->rowCopy()) && !(flags_ & `2`));
2007	}
2008	#ifndef CLP_ALL_ONE_FILE
2009	// These have to match ClpPrimalColumnSteepest version
2010	#define reference(i) (((reference[i>>5]>>(i&31))&1)!=0)
2011	#endif
2012	// Updates two arrays for steepest
2013	void
2014	ClpPackedMatrix::transposeTimes2(const ClpSimplex * model,
2015	const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
2016	const CoinIndexedVector * pi2,
2017	CoinIndexedVector * spare,
2018	double referenceIn, double devex,
2019	// Array for exact devex to say what is in reference framework
2020	unsigned int * reference,
2021	double * weights, double scaleFactor)
2022	{
2023	// put row of tableau in dj1
2024	double * pi = pi1->denseVector();
2025	int numberNonZero = `0`;
2026	int * index = dj1->getIndices();
2027	double * array = dj1->denseVector();
2028	int numberInRowArray = pi1->getNumElements();
2029	double zeroTolerance = model->zeroTolerance();
2030	bool packed = pi1->packedMode();
2031	// do by column
2032	int iColumn;
2033	// get matrix data pointers
2034	const int * row = matrix_->getIndices();
2035	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2036	const double * elementByColumn = matrix_->getElements();
2037	const double * rowScale = model->rowScale();
2038	assert (!spare->getNumElements());
2039	assert (numberActiveColumns_ > `0`);
2040	double * piWeight = pi2->denseVector();
2041	assert (!pi2->packedMode());
2042	bool killDjs = (scaleFactor == `0.0`);
2043	if (!scaleFactor)
2044	scaleFactor = `1.0`;
2045	if (packed) {
2046	// need to expand pi into y
2047	assert(spare->capacity() >= model->numberRows());
2048	double * piOld = pi;
2049	pi = spare->denseVector();
2050	const int * whichRow = pi1->getIndices();
2051	int i;
2052	ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
2053	if (rowScale && scaledMatrix) {
2054	rowScale = NULL;
2055	// get matrix data pointers
2056	row = scaledMatrix->getIndices();
2057	columnStart = scaledMatrix->getVectorStarts();
2058	elementByColumn = scaledMatrix->getElements();
2059	}
2060	if (!rowScale) {
2061	// modify pi so can collapse to one loop
2062	for (i = `0`; i < numberInRowArray; i++) {
2063	int iRow = whichRow[i];
2064	pi[iRow] = piOld[i];
2065	}
2066	if (!columnCopy_) {
2067	CoinBigIndex j;
2068	CoinBigIndex end = columnStart[`0`];
2069	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
2070	CoinBigIndex start = end;
2071	end = columnStart[iColumn+`1`];
2072	ClpSimplex::Status status = model->getStatus(iColumn);
2073	if (status == ClpSimplex::basic \|\| status == ClpSimplex::isFixed) continue;
2074	double value = `0.0`;
2075	for (j = start; j < end; j++) {
2076	int iRow = row[j];
2077	value -= pi[iRow] * elementByColumn[j];
2078	}
2079	if (fabs(value) > zeroTolerance) {
2080	// and do other array
2081	double modification = `0.0`;
2082	for (j = start; j < end; j++) {
2083	int iRow = row[j];
2084	modification += piWeight[iRow] * elementByColumn[j];
2085	}
2086	double thisWeight = weights[iColumn];
2087	double pivot = value * scaleFactor;
2088	double pivotSquared = pivot * pivot;
2089	thisWeight += pivotSquared * devex + pivot * modification;
2090	if (thisWeight < DEVEX_TRY_NORM) {
2091	if (referenceIn < `0.0`) {
2092	// steepest
2093	thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2094	} else {
2095	// exact
2096	thisWeight = referenceIn * pivotSquared;
2097	if (reference(iColumn))
2098	thisWeight += `1.0`;
2099	thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2100	}
2101	}
2102	weights[iColumn] = thisWeight;
2103	if (!killDjs) {
2104	array[numberNonZero] = value;
2105	index[numberNonZero++] = iColumn;
2106	}
2107	}
2108	}
2109	} else {
2110	// use special column copy
2111	// reset back
2112	if (killDjs)
2113	scaleFactor = `0.0`;
2114	columnCopy_->transposeTimes2(model, pi, dj1, piWeight, referenceIn, devex,
2115	reference, weights, scaleFactor);
2116	numberNonZero = dj1->getNumElements();
2117	}
2118	} else {
2119	// scaled
2120	// modify pi so can collapse to one loop
2121	for (i = `0`; i < numberInRowArray; i++) {
2122	int iRow = whichRow[i];
2123	pi[iRow] = piOld[i] * rowScale[iRow];
2124	}
2125	// can also scale piWeight as not used again
2126	int numberWeight = pi2->getNumElements();
2127	const int * indexWeight = pi2->getIndices();
2128	for (i = `0`; i < numberWeight; i++) {
2129	int iRow = indexWeight[i];
2130	piWeight[iRow] *= rowScale[iRow];
2131	}
2132	if (!columnCopy_) {
2133	const double * columnScale = model->columnScale();
2134	CoinBigIndex j;
2135	CoinBigIndex end = columnStart[`0`];
2136	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
2137	CoinBigIndex start = end;
2138	end = columnStart[iColumn+`1`];
2139	ClpSimplex::Status status = model->getStatus(iColumn);
2140	if (status == ClpSimplex::basic \|\| status == ClpSimplex::isFixed) continue;
2141	double scale = columnScale[iColumn];
2142	double value = `0.0`;
2143	for (j = start; j < end; j++) {
2144	int iRow = row[j];
2145	value -= pi[iRow] * elementByColumn[j];
2146	}
2147	value *= scale;
2148	if (fabs(value) > zeroTolerance) {
2149	double modification = `0.0`;
2150	for (j = start; j < end; j++) {
2151	int iRow = row[j];
2152	modification += piWeight[iRow] * elementByColumn[j];
2153	}
2154	modification *= scale;
2155	double thisWeight = weights[iColumn];
2156	double pivot = value * scaleFactor;
2157	double pivotSquared = pivot * pivot;
2158	thisWeight += pivotSquared * devex + pivot * modification;
2159	if (thisWeight < DEVEX_TRY_NORM) {
2160	if (referenceIn < `0.0`) {
2161	// steepest
2162	thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2163	} else {
2164	// exact
2165	thisWeight = referenceIn * pivotSquared;
2166	if (reference(iColumn))
2167	thisWeight += `1.0`;
2168	thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2169	}
2170	}
2171	weights[iColumn] = thisWeight;
2172	if (!killDjs) {
2173	array[numberNonZero] = value;
2174	index[numberNonZero++] = iColumn;
2175	}
2176	}
2177	}
2178	} else {
2179	// use special column copy
2180	// reset back
2181	if (killDjs)
2182	scaleFactor = `0.0`;
2183	columnCopy_->transposeTimes2(model, pi, dj1, piWeight, referenceIn, devex,
2184	reference, weights, scaleFactor);
2185	numberNonZero = dj1->getNumElements();
2186	}
2187	}
2188	// zero out
2189	int numberRows = model->numberRows();
2190	if (numberInRowArray * `4` < numberRows) {
2191	for (i = `0`; i < numberInRowArray; i++) {
2192	int iRow = whichRow[i];
2193	pi[iRow] = `0.0`;
2194	}
2195	} else {
2196	CoinZeroN(pi, numberRows);
2197	}
2198	} else {
2199	if (!rowScale) {
2200	CoinBigIndex j;
2201	CoinBigIndex end = columnStart[`0`];
2202	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
2203	CoinBigIndex start = end;
2204	end = columnStart[iColumn+`1`];
2205	ClpSimplex::Status status = model->getStatus(iColumn);
2206	if (status == ClpSimplex::basic \|\| status == ClpSimplex::isFixed) continue;
2207	double value = `0.0`;
2208	for (j = start; j < end; j++) {
2209	int iRow = row[j];
2210	value -= pi[iRow] * elementByColumn[j];
2211	}
2212	if (fabs(value) > zeroTolerance) {
2213	// and do other array
2214	double modification = `0.0`;
2215	for (j = start; j < end; j++) {
2216	int iRow = row[j];
2217	modification += piWeight[iRow] * elementByColumn[j];
2218	}
2219	double thisWeight = weights[iColumn];
2220	double pivot = value * scaleFactor;
2221	double pivotSquared = pivot * pivot;
2222	thisWeight += pivotSquared * devex + pivot * modification;
2223	if (thisWeight < DEVEX_TRY_NORM) {
2224	if (referenceIn < `0.0`) {
2225	// steepest
2226	thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2227	} else {
2228	// exact
2229	thisWeight = referenceIn * pivotSquared;
2230	if (reference(iColumn))
2231	thisWeight += `1.0`;
2232	thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2233	}
2234	}
2235	weights[iColumn] = thisWeight;
2236	if (!killDjs) {
2237	array[iColumn] = value;
2238	index[numberNonZero++] = iColumn;
2239	}
2240	}
2241	}
2242	} else {
2243	#ifdef CLP_INVESTIGATE
2244	if (model->clpScaledMatrix())
2245	printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
2246	#endif
2247	// scaled
2248	// can also scale piWeight as not used again
2249	int numberWeight = pi2->getNumElements();
2250	const int * indexWeight = pi2->getIndices();
2251	for (int i = `0`; i < numberWeight; i++) {
2252	int iRow = indexWeight[i];
2253	piWeight[iRow] *= rowScale[iRow];
2254	}
2255	const double * columnScale = model->columnScale();
2256	CoinBigIndex j;
2257	CoinBigIndex end = columnStart[`0`];
2258	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
2259	CoinBigIndex start = end;
2260	end = columnStart[iColumn+`1`];
2261	ClpSimplex::Status status = model->getStatus(iColumn);
2262	if (status == ClpSimplex::basic \|\| status == ClpSimplex::isFixed) continue;
2263	double scale = columnScale[iColumn];
2264	double value = `0.0`;
2265	for (j = start; j < end; j++) {
2266	int iRow = row[j];
2267	value -= pi[iRow] * elementByColumn[j] * rowScale[iRow];
2268	}
2269	value *= scale;
2270	if (fabs(value) > zeroTolerance) {
2271	double modification = `0.0`;
2272	for (j = start; j < end; j++) {
2273	int iRow = row[j];
2274	modification += piWeight[iRow] * elementByColumn[j];
2275	}
2276	modification *= scale;
2277	double thisWeight = weights[iColumn];
2278	double pivot = value * scaleFactor;
2279	double pivotSquared = pivot * pivot;
2280	thisWeight += pivotSquared * devex + pivot * modification;
2281	if (thisWeight < DEVEX_TRY_NORM) {
2282	if (referenceIn < `0.0`) {
2283	// steepest
2284	thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2285	} else {
2286	// exact
2287	thisWeight = referenceIn * pivotSquared;
2288	if (reference(iColumn))
2289	thisWeight += `1.0`;
2290	thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2291	}
2292	}
2293	weights[iColumn] = thisWeight;
2294	if (!killDjs) {
2295	array[iColumn] = value;
2296	index[numberNonZero++] = iColumn;
2297	}
2298	}
2299	}
2300	}
2301	}
2302	dj1->setNumElements(numberNonZero);
2303	spare->setNumElements(`0`);
2304	if (packed)
2305	dj1->setPackedMode(true);
2306	}
2307	// Updates second array for steepest and does devex weights
2308	void
2309	ClpPackedMatrix::subsetTimes2(const ClpSimplex * model,
2310	CoinIndexedVector * dj1,
2311	const CoinIndexedVector * pi2, CoinIndexedVector *,
2312	double referenceIn, double devex,
2313	// Array for exact devex to say what is in reference framework
2314	unsigned int * reference,
2315	double * weights, double scaleFactor)
2316	{
2317	int number = dj1->getNumElements();
2318	const int * index = dj1->getIndices();
2319	double * array = dj1->denseVector();
2320	assert( dj1->packedMode());
2321
2322	// get matrix data pointers
2323	const int * row = matrix_->getIndices();
2324	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2325	const int * columnLength = matrix_->getVectorLengths();
2326	const double * elementByColumn = matrix_->getElements();
2327	const double * rowScale = model->rowScale();
2328	double * piWeight = pi2->denseVector();
2329	bool killDjs = (scaleFactor == `0.0`);
2330	if (!scaleFactor)
2331	scaleFactor = `1.0`;
2332	if (!rowScale) {
2333	for (int k = `0`; k < number; k++) {
2334	int iColumn = index[k];
2335	double pivot = array[k] * scaleFactor;
2336	if (killDjs)
2337	array[k] = `0.0`;
2338	// and do other array
2339	double modification = `0.0`;
2340	for (CoinBigIndex j = columnStart[iColumn];
2341	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2342	int iRow = row[j];
2343	modification += piWeight[iRow] * elementByColumn[j];
2344	}
2345	double thisWeight = weights[iColumn];
2346	double pivotSquared = pivot * pivot;
2347	thisWeight += pivotSquared * devex + pivot * modification;
2348	if (thisWeight < DEVEX_TRY_NORM) {
2349	if (referenceIn < `0.0`) {
2350	// steepest
2351	thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2352	} else {
2353	// exact
2354	thisWeight = referenceIn * pivotSquared;
2355	if (reference(iColumn))
2356	thisWeight += `1.0`;
2357	thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2358	}
2359	}
2360	weights[iColumn] = thisWeight;
2361	}
2362	} else {
2363	#ifdef CLP_INVESTIGATE
2364	if (model->clpScaledMatrix())
2365	printf("scaledMatrix_ at %d of ClpPackedMatrix\n", __LINE__);
2366	#endif
2367	// scaled
2368	const double * columnScale = model->columnScale();
2369	for (int k = `0`; k < number; k++) {
2370	int iColumn = index[k];
2371	double pivot = array[k] * scaleFactor;
2372	double scale = columnScale[iColumn];
2373	if (killDjs)
2374	array[k] = `0.0`;
2375	// and do other array
2376	double modification = `0.0`;
2377	for (CoinBigIndex j = columnStart[iColumn];
2378	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2379	int iRow = row[j];
2380	modification += piWeight[iRow] * elementByColumn[j] * rowScale[iRow];
2381	}
2382	double thisWeight = weights[iColumn];
2383	modification *= scale;
2384	double pivotSquared = pivot * pivot;
2385	thisWeight += pivotSquared * devex + pivot * modification;
2386	if (thisWeight < DEVEX_TRY_NORM) {
2387	if (referenceIn < `0.0`) {
2388	// steepest
2389	thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2390	} else {
2391	// exact
2392	thisWeight = referenceIn * pivotSquared;
2393	if (reference(iColumn))
2394	thisWeight += `1.0`;
2395	thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2396	}
2397	}
2398	weights[iColumn] = thisWeight;
2399	}
2400	}
2401	}
2402	/// returns number of elements in column part of basis,
2403	CoinBigIndex
2404	ClpPackedMatrix::countBasis( const int * whichColumn,
2405	int & numberColumnBasic)
2406	{
2407	const int * columnLength = matrix_->getVectorLengths();
2408	int i;
2409	CoinBigIndex numberElements = `0`;
2410	// just count - can be over so ignore zero problem
2411	for (i = `0`; i < numberColumnBasic; i++) {
2412	int iColumn = whichColumn[i];
2413	numberElements += columnLength[iColumn];
2414	}
2415	return numberElements;
2416	}
2417	void
2418	ClpPackedMatrix::fillBasis(ClpSimplex * model,
2419	const int * COIN_RESTRICT whichColumn,
2420	int & numberColumnBasic,
2421	int * COIN_RESTRICT indexRowU,
2422	int * COIN_RESTRICT start,
2423	int * COIN_RESTRICT rowCount,
2424	int * COIN_RESTRICT columnCount,
2425	CoinFactorizationDouble * COIN_RESTRICT elementU)
2426	{
2427	const int * COIN_RESTRICT columnLength = matrix_->getVectorLengths();
2428	int i;
2429	CoinBigIndex numberElements = start[`0`];
2430	// fill
2431	const CoinBigIndex * COIN_RESTRICT columnStart = matrix_->getVectorStarts();
2432	const double * COIN_RESTRICT rowScale = model->rowScale();
2433	const int * COIN_RESTRICT row = matrix_->getIndices();
2434	const double * COIN_RESTRICT elementByColumn = matrix_->getElements();
2435	ClpPackedMatrix * scaledMatrix = model->clpScaledMatrix();
2436	if (scaledMatrix && true) {
2437	columnLength = scaledMatrix->matrix_->getVectorLengths();
2438	columnStart = scaledMatrix->matrix_->getVectorStarts();
2439	rowScale = NULL;
2440	row = scaledMatrix->matrix_->getIndices();
2441	elementByColumn = scaledMatrix->matrix_->getElements();
2442	}
2443	if ((flags_ & `1`) == `0`) {
2444	if (!rowScale) {
2445	// no scaling
2446	for (i = `0`; i < numberColumnBasic; i++) {
2447	int iColumn = whichColumn[i];
2448	int length = columnLength[iColumn];
2449	CoinBigIndex startThis = columnStart[iColumn];
2450	columnCount[i] = length;
2451	CoinBigIndex endThis = startThis + length;
2452	for (CoinBigIndex j = startThis; j < endThis; j++) {
2453	int iRow = row[j];
2454	indexRowU[numberElements] = iRow;
2455	rowCount[iRow]++;
2456	assert (elementByColumn[j]);
2457	elementU[numberElements++] = elementByColumn[j];
2458	}
2459	start[i+`1`] = numberElements;
2460	}
2461	} else {
2462	// scaling
2463	const double * COIN_RESTRICT columnScale = model->columnScale();
2464	for (i = `0`; i < numberColumnBasic; i++) {
2465	int iColumn = whichColumn[i];
2466	double scale = columnScale[iColumn];
2467	int length = columnLength[iColumn];
2468	CoinBigIndex startThis = columnStart[iColumn];
2469	columnCount[i] = length;
2470	CoinBigIndex endThis = startThis + length;
2471	for (CoinBigIndex j = startThis; j < endThis; j++) {
2472	int iRow = row[j];
2473	indexRowU[numberElements] = iRow;
2474	rowCount[iRow]++;
2475	assert (elementByColumn[j]);
2476	elementU[numberElements++] =
2477	elementByColumn[j] * scale * rowScale[iRow];
2478	}
2479	start[i+`1`] = numberElements;
2480	}
2481	}
2482	} else {
2483	// there are zero elements so need to look more closely
2484	if (!rowScale) {
2485	// no scaling
2486	for (i = `0`; i < numberColumnBasic; i++) {
2487	int iColumn = whichColumn[i];
2488	CoinBigIndex j;
2489	for (j = columnStart[iColumn];
2490	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2491	double value = elementByColumn[j];
2492	if (value) {
2493	int iRow = row[j];
2494	indexRowU[numberElements] = iRow;
2495	rowCount[iRow]++;
2496	elementU[numberElements++] = value;
2497	}
2498	}
2499	start[i+`1`] = numberElements;
2500	columnCount[i] = numberElements - start[i];
2501	}
2502	} else {
2503	// scaling
2504	const double * columnScale = model->columnScale();
2505	for (i = `0`; i < numberColumnBasic; i++) {
2506	int iColumn = whichColumn[i];
2507	CoinBigIndex j;
2508	double scale = columnScale[iColumn];
2509	for (j = columnStart[iColumn];
2510	j < columnStart[iColumn] + columnLength[i]; j++) {
2511	double value = elementByColumn[j];
2512	if (value) {
2513	int iRow = row[j];
2514	indexRowU[numberElements] = iRow;
2515	rowCount[iRow]++;
2516	elementU[numberElements++] = value * scale * rowScale[iRow];
2517	}
2518	}
2519	start[i+`1`] = numberElements;
2520	columnCount[i] = numberElements - start[i];
2521	}
2522	}
2523	}
2524	}
2525	#if 0
2526	int
2527	ClpPackedMatrix::scale2(ClpModel * model) const
2528	{
2529	ClpSimplex * baseModel = NULL;
2530	#ifndef NDEBUG
2531	//checkFlags();
2532	#endif
2533	int numberRows = model->numberRows();
2534	int numberColumns = matrix_->getNumCols();
2535	model->setClpScaledMatrix(NULL); // get rid of any scaled matrix
2536	// If empty - return as sanityCheck will trap
2537	if (!numberRows \|\| !numberColumns) {
2538	model->setRowScale(NULL);
2539	model->setColumnScale(NULL);
2540	return `1`;
2541	}
2542	ClpMatrixBase * rowCopyBase = model->rowCopy();
2543	double * rowScale;
2544	double * columnScale;
2545	//assert (!model->rowScale());
2546	bool arraysExist;
2547	double * inverseRowScale = NULL;
2548	double * inverseColumnScale = NULL;
2549	if (!model->rowScale()) {
2550	rowScale = new double [numberRows*`2`];
2551	columnScale = new double [numberColumns*`2`];
2552	inverseRowScale = rowScale + numberRows;
2553	inverseColumnScale = columnScale + numberColumns;
2554	arraysExist = false;
2555	} else {
2556	rowScale = model->mutableRowScale();
2557	columnScale = model->mutableColumnScale();
2558	inverseRowScale = model->mutableInverseRowScale();
2559	inverseColumnScale = model->mutableInverseColumnScale();
2560	arraysExist = true;
2561	}
2562	assert (inverseRowScale == rowScale + numberRows);
2563	assert (inverseColumnScale == columnScale + numberColumns);
2564	// we are going to mark bits we are interested in
2565	char * usefulRow = new char [numberRows];
2566	char * usefulColumn = new char [numberColumns];
2567	double * rowLower = model->rowLower();
2568	double * rowUpper = model->rowUpper();
2569	double * columnLower = model->columnLower();
2570	double * columnUpper = model->columnUpper();
2571	int iColumn, iRow;
2572	//#define LEAVE_FIXED
2573	// mark free rows
2574	for (iRow = `0`; iRow < numberRows; iRow++) {
2575	#if 0 //ndef LEAVE_FIXED
2576	if (rowUpper[iRow] < `1.0e20` \|\|
2577	rowLower[iRow] > -`1.0e20`)
2578	usefulRow[iRow] = `1`;
2579	else
2580	usefulRow[iRow] = `0`;
2581	#else
2582	usefulRow[iRow] = `1`;
2583	#endif
2584	}
2585	// mark empty and fixed columns
2586	// also see if worth scaling
2587	assert (model->scalingFlag() <= `4`);
2588	// scale_stats[model->scalingFlag()]++;
2589	double largest = `0.0`;
2590	double smallest = `1.0e50`;
2591	// get matrix data pointers
2592	int * row = matrix_->getMutableIndices();
2593	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
2594	int * columnLength = matrix_->getMutableVectorLengths();
2595	double * elementByColumn = matrix_->getMutableElements();
2596	bool deletedElements = false;
2597	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
2598	CoinBigIndex j;
2599	char useful = `0`;
2600	bool deleteSome = false;
2601	int start = columnStart[iColumn];
2602	int end = start + columnLength[iColumn];
2603	#ifndef LEAVE_FIXED
2604	if (columnUpper[iColumn] >
2605	columnLower[iColumn] + `1.0e-12`) {
2606	#endif
2607	for (j = start; j < end; j++) {
2608	iRow = row[j];
2609	double value = fabs(elementByColumn[j]);
2610	if (value > `1.0e-20`) {
2611	if(usefulRow[iRow]) {
2612	useful = `1`;
2613	largest = CoinMax(largest, fabs(elementByColumn[j]));
2614	smallest = CoinMin(smallest, fabs(elementByColumn[j]));
2615	}
2616	} else {
2617	// small
2618	deleteSome = true;
2619	}
2620	}
2621	#ifndef LEAVE_FIXED
2622	} else {
2623	// just check values
2624	for (j = start; j < end; j++) {
2625	double value = fabs(elementByColumn[j]);
2626	if (value <= `1.0e-20`) {
2627	// small
2628	deleteSome = true;
2629	}
2630	}
2631	}
2632	#endif
2633	usefulColumn[iColumn] = useful;
2634	if (deleteSome) {
2635	deletedElements = true;
2636	CoinBigIndex put = start;
2637	for (j = start; j < end; j++) {
2638	double value = elementByColumn[j];
2639	if (fabs(value) > `1.0e-20`) {
2640	row[put] = row[j];
2641	elementByColumn[put++] = value;
2642	}
2643	}
2644	columnLength[iColumn] = put - start;
2645	}
2646	}
2647	model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL, *model->messagesPointer())
2648	<< smallest << largest
2649	<< CoinMessageEol;
2650	if (smallest >= `0.5` && largest <= `2.0` && !deletedElements) {
2651	// don't bother scaling
2652	model->messageHandler()->message(CLP_PACKEDSCALE_FORGET, *model->messagesPointer())
2653	<< CoinMessageEol;
2654	if (!arraysExist) {
2655	delete [] rowScale;
2656	delete [] columnScale;
2657	} else {
2658	model->setRowScale(NULL);
2659	model->setColumnScale(NULL);
2660	}
2661	delete [] usefulRow;
2662	delete [] usefulColumn;
2663	return `1`;
2664	} else {
2665	#ifdef CLP_INVESTIGATE
2666	if (deletedElements)
2667	printf("DEL_ELS\n");
2668	#endif
2669	if (!rowCopyBase) {
2670	// temporary copy
2671	rowCopyBase = reverseOrderedCopy();
2672	} else if (deletedElements) {
2673	rowCopyBase = reverseOrderedCopy();
2674	model->setNewRowCopy(rowCopyBase);
2675	}
2676	#ifndef NDEBUG
2677	ClpPackedMatrix* rowCopy =
2678	dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
2679	// Make sure it is really a ClpPackedMatrix
2680	assert (rowCopy != NULL);
2681	#else
2682	ClpPackedMatrix* rowCopy =
2683	static_cast< ClpPackedMatrix*>(rowCopyBase);
2684	#endif
2685
2686	const int * column = rowCopy->getIndices();
2687	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
2688	const double * element = rowCopy->getElements();
2689	// need to scale
2690	if (largest > `1.0e13` * smallest) {
2691	// safer to have smaller zero tolerance
2692	double ratio = smallest / largest;
2693	ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
2694	double newTolerance = CoinMax(ratio * `0.5`, `1.0e-18`);
2695	if (simplex->zeroTolerance() > newTolerance)
2696	simplex->setZeroTolerance(newTolerance);
2697	}
2698	int scalingMethod = model->scalingFlag();
2699	if (scalingMethod == `4`) {
2700	// As auto
2701	scalingMethod = `3`;
2702	}
2703	// and see if there any empty rows
2704	for (iRow = `0`; iRow < numberRows; iRow++) {
2705	if (usefulRow[iRow]) {
2706	CoinBigIndex j;
2707	int useful = `0`;
2708	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
2709	int iColumn = column[j];
2710	if (usefulColumn[iColumn]) {
2711	useful = `1`;
2712	break;
2713	}
2714	}
2715	usefulRow[iRow] = static_cast<char>(useful);
2716	}
2717	}
2718	double savedOverallRatio = `0.0`;
2719	double tolerance = `5.0` * model->primalTolerance();
2720	double overallLargest = -`1.0e-20`;
2721	double overallSmallest = `1.0e20`;
2722	bool finished = false;
2723	// if scalingMethod 3 then may change
2724	bool extraDetails = (model->logLevel() > `2`);
2725	while (!finished) {
2726	int numberPass = `3`;
2727	overallLargest = -`1.0e-20`;
2728	overallSmallest = `1.0e20`;
2729	if (!baseModel) {
2730	ClpFillN ( rowScale, numberRows, `1.0`);
2731	ClpFillN ( columnScale, numberColumns, `1.0`);
2732	} else {
2733	// Copy scales and do quick scale on extra rows
2734	// Then just one? pass
2735	assert (numberColumns == baseModel->numberColumns());
2736	int numberRows2 = baseModel->numberRows();
2737	assert (numberRows >= numberRows2);
2738	assert (baseModel->rowScale());
2739	CoinMemcpyN(baseModel->rowScale(), numberRows2, rowScale);
2740	CoinMemcpyN(baseModel->columnScale(), numberColumns, columnScale);
2741	if (numberRows > numberRows2) {
2742	numberPass = `1`;
2743	// do some scaling
2744	if (scalingMethod == `1` \|\| scalingMethod == `3`) {
2745	// Maximum in each row
2746	for (iRow = numberRows2; iRow < numberRows; iRow++) {
2747	if (usefulRow[iRow]) {
2748	CoinBigIndex j;
2749	largest = `1.0e-10`;
2750	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
2751	int iColumn = column[j];
2752	if (usefulColumn[iColumn]) {
2753	double value = fabs(element[j] * columnScale[iColumn]);
2754	largest = CoinMax(largest, value);
2755	assert (largest < `1.0e40`);
2756	}
2757	}
2758	rowScale[iRow] = `1.0` / largest;
2759	#ifdef COIN_DEVELOP
2760	if (extraDetails) {
2761	overallLargest = CoinMax(overallLargest, largest);
2762	overallSmallest = CoinMin(overallSmallest, largest);
2763	}
2764	#endif
2765	}
2766	}
2767	} else {
2768	overallLargest = `0.0`;
2769	overallSmallest = `1.0e50`;
2770	// Geometric mean on row scales
2771	for (iRow = `0`; iRow < numberRows; iRow++) {
2772	if (usefulRow[iRow]) {
2773	CoinBigIndex j;
2774	largest = `1.0e-20`;
2775	smallest = `1.0e50`;
2776	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
2777	int iColumn = column[j];
2778	if (usefulColumn[iColumn]) {
2779	double value = fabs(element[j]);
2780	value *= columnScale[iColumn];
2781	largest = CoinMax(largest, value);
2782	smallest = CoinMin(smallest, value);
2783	}
2784	}
2785	if (iRow >= numberRows2) {
2786	rowScale[iRow] = `1.0` / sqrt(smallest * largest);
2787	//rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
2788	}
2789	#ifdef COIN_DEVELOP
2790	if (extraDetails) {
2791	overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
2792	overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
2793	}
2794	#endif
2795	}
2796	}
2797	}
2798	} else {
2799	// just use
2800	numberPass = `0`;
2801	}
2802	}
2803	if (!baseModel && (scalingMethod == `1` \|\| scalingMethod == `3`)) {
2804	// Maximum in each row
2805	for (iRow = `0`; iRow < numberRows; iRow++) {
2806	if (usefulRow[iRow]) {
2807	CoinBigIndex j;
2808	largest = `1.0e-10`;
2809	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
2810	int iColumn = column[j];
2811	if (usefulColumn[iColumn]) {
2812	double value = fabs(element[j]);
2813	largest = CoinMax(largest, value);
2814	assert (largest < `1.0e40`);
2815	}
2816	}
2817	rowScale[iRow] = `1.0` / largest;
2818	#ifdef COIN_DEVELOP
2819	if (extraDetails) {
2820	overallLargest = CoinMax(overallLargest, largest);
2821	overallSmallest = CoinMin(overallSmallest, largest);
2822	}
2823	#endif
2824	}
2825	}
2826	} else {
2827	#ifdef USE_OBJECTIVE
2828	// This will be used to help get scale factors
2829	double * objective = new double[numberColumns];
2830	CoinMemcpyN(model->costRegion(`1`), numberColumns, objective);
2831	double objScale = `1.0`;
2832	#endif
2833	while (numberPass) {
2834	overallLargest = `0.0`;
2835	overallSmallest = `1.0e50`;
2836	numberPass--;
2837	// Geometric mean on row scales
2838	for (iRow = `0`; iRow < numberRows; iRow++) {
2839	if (usefulRow[iRow]) {
2840	CoinBigIndex j;
2841	largest = `1.0e-20`;
2842	smallest = `1.0e50`;
2843	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
2844	int iColumn = column[j];
2845	if (usefulColumn[iColumn]) {
2846	double value = fabs(element[j]);
2847	value *= columnScale[iColumn];
2848	largest = CoinMax(largest, value);
2849	smallest = CoinMin(smallest, value);
2850	}
2851	}
2852
2853	rowScale[iRow] = `1.0` / sqrt(smallest * largest);
2854	//rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
2855	if (extraDetails) {
2856	overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
2857	overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
2858	}
2859	}
2860	}
2861	#ifdef USE_OBJECTIVE
2862	largest = `1.0e-20`;
2863	smallest = `1.0e50`;
2864	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
2865	if (usefulColumn[iColumn]) {
2866	double value = fabs(objective[iColumn]);
2867	value *= columnScale[iColumn];
2868	largest = CoinMax(largest, value);
2869	smallest = CoinMin(smallest, value);
2870	}
2871	}
2872	objScale = `1.0` / sqrt(smallest * largest);
2873	#endif
2874	model->messageHandler()->message(CLP_PACKEDSCALE_WHILE, *model->messagesPointer())
2875	<< overallSmallest
2876	<< overallLargest
2877	<< CoinMessageEol;
2878	// skip last column round
2879	if (numberPass == `1`)
2880	break;
2881	// Geometric mean on column scales
2882	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
2883	if (usefulColumn[iColumn]) {
2884	CoinBigIndex j;
2885	largest = `1.0e-20`;
2886	smallest = `1.0e50`;
2887	for (j = columnStart[iColumn];
2888	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2889	iRow = row[j];
2890	double value = fabs(elementByColumn[j]);
2891	if (usefulRow[iRow]) {
2892	value *= rowScale[iRow];
2893	largest = CoinMax(largest, value);
2894	smallest = CoinMin(smallest, value);
2895	}
2896	}
2897	#ifdef USE_OBJECTIVE
2898	if (fabs(objective[iColumn]) > `1.0e-20`) {
2899	double value = fabs(objective[iColumn]) * objScale;
2900	largest = CoinMax(largest, value);
2901	smallest = CoinMin(smallest, value);
2902	}
2903	#endif
2904	columnScale[iColumn] = `1.0` / sqrt(smallest * largest);
2905	//columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
2906	}
2907	}
2908	}
2909	#ifdef USE_OBJECTIVE
2910	delete [] objective;
2911	printf("obj scale %g - use it if you want\n", objScale);
2912	#endif
2913	}
2914	// If ranges will make horrid then scale
2915	for (iRow = `0`; iRow < numberRows; iRow++) {
2916	if (usefulRow[iRow]) {
2917	double difference = rowUpper[iRow] - rowLower[iRow];
2918	double scaledDifference = difference * rowScale[iRow];
2919	if (scaledDifference > tolerance && scaledDifference < `1.0e-4`) {
2920	// make gap larger
2921	rowScale[iRow] *= `1.0e-4` / scaledDifference;
2922	rowScale[iRow] = CoinMax(`1.0e-10`, CoinMin(`1.0e10`, rowScale[iRow]));
2923	//printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
2924	// scaledDifference,differencerowScale[iRow]);*
2925	}
2926	}
2927	}
2928	// final pass to scale columns so largest is reasonable
2929	// See what smallest will be if largest is 1.0
2930	overallSmallest = `1.0e50`;
2931	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
2932	if (usefulColumn[iColumn]) {
2933	CoinBigIndex j;
2934	largest = `1.0e-20`;
2935	smallest = `1.0e50`;
2936	for (j = columnStart[iColumn];
2937	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
2938	iRow = row[j];
2939	if(elementByColumn[j] && usefulRow[iRow]) {
2940	double value = fabs(elementByColumn[j] * rowScale[iRow]);
2941	largest = CoinMax(largest, value);
2942	smallest = CoinMin(smallest, value);
2943	}
2944	}
2945	if (overallSmallest * largest > smallest)
2946	overallSmallest = smallest / largest;
2947	}
2948	}
2949	if (scalingMethod == `1` \|\| scalingMethod == `2`) {
2950	finished = true;
2951	} else if (savedOverallRatio == `0.0` && scalingMethod != `4`) {
2952	savedOverallRatio = overallSmallest;
2953	scalingMethod = `4`;
2954	} else {
2955	assert (scalingMethod == `4`);
2956	if (overallSmallest > `2.0` * savedOverallRatio) {
2957	finished = true; // geometric was better
2958	if (model->scalingFlag() == `4`) {
2959	// If in Branch and bound change
2960	if ((model->specialOptions() & `1024`) != `0`) {
2961	model->scaling(`2`);
2962	}
2963	}
2964	} else {
2965	scalingMethod = `1`; // redo equilibrium
2966	if (model->scalingFlag() == `4`) {
2967	// If in Branch and bound change
2968	if ((model->specialOptions() & `1024`) != `0`) {
2969	model->scaling(`1`);
2970	}
2971	}
2972	}
2973	#if 0
2974	if (extraDetails) {
2975	if (finished)
2976	printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
2977	savedOverallRatio, overallSmallest);
2978	else
2979	printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
2980	savedOverallRatio, overallSmallest);
2981	}
2982	#endif
2983	}
2984	}
2985	//#define RANDOMIZE
2986	#ifdef RANDOMIZE
2987	// randomize by up to 10%
2988	for (iRow = `0`; iRow < numberRows; iRow++) {
2989	double value = `0.5` - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
2990	rowScale[iRow] = (`1.0` + `0.1` value);
2991	}
2992	#endif
2993	overallLargest = `1.0`;
2994	if (overallSmallest < `1.0e-1`)
2995	overallLargest = `1.0` / sqrt(overallSmallest);
2996	overallLargest = CoinMin(`100.0`, overallLargest);
2997	overallSmallest = `1.0e50`;
2998	//printf("scaling %d\n",model->scalingFlag());
2999	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
3000	if (columnUpper[iColumn] >
3001	columnLower[iColumn] + `1.0e-12`) {
3002	//if (usefulColumn[iColumn]) {
3003	CoinBigIndex j;
3004	largest = `1.0e-20`;
3005	smallest = `1.0e50`;
3006	for (j = columnStart[iColumn];
3007	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3008	iRow = row[j];
3009	if(elementByColumn[j] && usefulRow[iRow]) {
3010	double value = fabs(elementByColumn[j] * rowScale[iRow]);
3011	largest = CoinMax(largest, value);
3012	smallest = CoinMin(smallest, value);
3013	}
3014	}
3015	columnScale[iColumn] = overallLargest / largest;
3016	//columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
3017	#ifdef RANDOMIZE
3018	double value = `0.5` - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
3019	columnScale[iColumn] = (`1.0` + `0.1` value);
3020	#endif
3021	double difference = columnUpper[iColumn] - columnLower[iColumn];
3022	if (difference < `1.0e-5` * columnScale[iColumn]) {
3023	// make gap larger
3024	columnScale[iColumn] = difference / `1.0e-5`;
3025	//printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
3026	// scaledDifference,differencecolumnScale[iColumn]);*
3027	}
3028	double value = smallest * columnScale[iColumn];
3029	if (overallSmallest > value)
3030	overallSmallest = value;
3031	//overallSmallest = CoinMin(overallSmallest,smallestcolumnScale[iColumn]);*
3032	}
3033	}
3034	model->messageHandler()->message(CLP_PACKEDSCALE_FINAL, *model->messagesPointer())
3035	<< overallSmallest
3036	<< overallLargest
3037	<< CoinMessageEol;
3038	if (overallSmallest < `1.0e-13`) {
3039	// Change factorization zero tolerance
3040	double newTolerance = CoinMax(`1.0e-15` * (overallSmallest / `1.0e-13`),
3041	`1.0e-18`);
3042	ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
3043	if (simplex->factorization()->zeroTolerance() > newTolerance)
3044	simplex->factorization()->zeroTolerance(newTolerance);
3045	newTolerance = CoinMax(overallSmallest * `0.5`, `1.0e-18`);
3046	simplex->setZeroTolerance(newTolerance);
3047	}
3048	delete [] usefulRow;
3049	delete [] usefulColumn;
3050	#ifndef SLIM_CLP
3051	// If quadratic then make symmetric
3052	ClpObjective * obj = model->objectiveAsObject();
3053	#ifndef NO_RTTI
3054	ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
3055	#else
3056	ClpQuadraticObjective * quadraticObj = NULL;
3057	if (obj->type() == `2`)
3058	quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
3059	#endif
3060	if (quadraticObj) {
3061	if (!rowCopyBase) {
3062	// temporary copy
3063	rowCopyBase = reverseOrderedCopy();
3064	}
3065	#ifndef NDEBUG
3066	ClpPackedMatrix* rowCopy =
3067	dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3068	// Make sure it is really a ClpPackedMatrix
3069	assert (rowCopy != NULL);
3070	#else
3071	ClpPackedMatrix* rowCopy =
3072	static_cast< ClpPackedMatrix*>(rowCopyBase);
3073	#endif
3074	const int * column = rowCopy->getIndices();
3075	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3076	CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
3077	int numberXColumns = quadratic->getNumCols();
3078	if (numberXColumns < numberColumns) {
3079	// we assume symmetric
3080	int numberQuadraticColumns = `0`;
3081	int i;
3082	//const int columnQuadratic = quadratic->getIndices();*
3083	//const int columnQuadraticStart = quadratic->getVectorStarts();*
3084	const int * columnQuadraticLength = quadratic->getVectorLengths();
3085	for (i = `0`; i < numberXColumns; i++) {
3086	int length = columnQuadraticLength[i];
3087	#ifndef CORRECT_COLUMN_COUNTS
3088	length = `1`;
3089	#endif
3090	if (length)
3091	numberQuadraticColumns++;
3092	}
3093	int numberXRows = numberRows - numberQuadraticColumns;
3094	numberQuadraticColumns = `0`;
3095	for (i = `0`; i < numberXColumns; i++) {
3096	int length = columnQuadraticLength[i];
3097	#ifndef CORRECT_COLUMN_COUNTS
3098	length = `1`;
3099	#endif
3100	if (length) {
3101	rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
3102	numberQuadraticColumns++;
3103	}
3104	}
3105	int numberQuadraticRows = `0`;
3106	for (i = `0`; i < numberXRows; i++) {
3107	// See if any in row quadratic
3108	CoinBigIndex j;
3109	int numberQ = `0`;
3110	for (j = rowStart[i]; j < rowStart[i+`1`]; j++) {
3111	int iColumn = column[j];
3112	if (columnQuadraticLength[iColumn])
3113	numberQ++;
3114	}
3115	#ifndef CORRECT_ROW_COUNTS
3116	numberQ = `1`;
3117	#endif
3118	if (numberQ) {
3119	columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
3120	numberQuadraticRows++;
3121	}
3122	}
3123	// and make sure Sj okay
3124	for (iColumn = numberQuadraticRows + numberXColumns; iColumn < numberColumns; iColumn++) {
3125	CoinBigIndex j = columnStart[iColumn];
3126	assert(columnLength[iColumn] == `1`);
3127	int iRow = row[j];
3128	double value = fabs(elementByColumn[j] * rowScale[iRow]);
3129	columnScale[iColumn] = `1.0` / value;
3130	}
3131	}
3132	}
3133	#endif
3134	// make copy (could do faster by using previous values)
3135	// could just do partial
3136	for (iRow = `0`; iRow < numberRows; iRow++)
3137	inverseRowScale[iRow] = `1.0` / rowScale[iRow] ;
3138	for (iColumn = `0`; iColumn < numberColumns; iColumn++)
3139	inverseColumnScale[iColumn] = `1.0` / columnScale[iColumn] ;
3140	if (!arraysExist) {
3141	model->setRowScale(rowScale);
3142	model->setColumnScale(columnScale);
3143	}
3144	if (model->rowCopy()) {
3145	// need to replace row by row
3146	ClpPackedMatrix* rowCopy =
3147	static_cast< ClpPackedMatrix*>(model->rowCopy());
3148	double * element = rowCopy->getMutableElements();
3149	const int * column = rowCopy->getIndices();
3150	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3151	// scale row copy
3152	for (iRow = `0`; iRow < numberRows; iRow++) {
3153	CoinBigIndex j;
3154	double scale = rowScale[iRow];
3155	double * elementsInThisRow = element + rowStart[iRow];
3156	const int * columnsInThisRow = column + rowStart[iRow];
3157	int number = rowStart[iRow+`1`] - rowStart[iRow];
3158	assert (number <= numberColumns);
3159	for (j = `0`; j < number; j++) {
3160	int iColumn = columnsInThisRow[j];
3161	elementsInThisRow[j] = scale columnScale[iColumn];
3162	}
3163	}
3164	if ((model->specialOptions() & `262144`) != `0`) {
3165	//if ((model->specialOptions()&(COIN_CBC_USING_CLP\|16384))!=0) {
3166	//if (model->inCbcBranchAndBound()&&false) {
3167	// copy without gaps
3168	CoinPackedMatrix * scaledMatrix = new CoinPackedMatrix(*matrix_, `0`, `0`);
3169	ClpPackedMatrix * scaled = new ClpPackedMatrix(scaledMatrix);
3170	model->setClpScaledMatrix(scaled);
3171	// get matrix data pointers
3172	const int * row = scaledMatrix->getIndices();
3173	const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
3174	#ifndef NDEBUG
3175	const int * columnLength = scaledMatrix->getVectorLengths();
3176	#endif
3177	double * elementByColumn = scaledMatrix->getMutableElements();
3178	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
3179	CoinBigIndex j;
3180	double scale = columnScale[iColumn];
3181	assert (columnStart[iColumn+`1`] == columnStart[iColumn] + columnLength[iColumn]);
3182	for (j = columnStart[iColumn];
3183	j < columnStart[iColumn+`1`]; j++) {
3184	int iRow = row[j];
3185	elementByColumn[j] = scale rowScale[iRow];
3186	}
3187	}
3188	} else {
3189	//printf("not in b&b\n");
3190	}
3191	} else {
3192	// no row copy
3193	delete rowCopyBase;
3194	}
3195	return `0`;
3196	}
3197	}
3198	#endif
3199	//#define SQRT_ARRAY
3200	#ifdef SQRT_ARRAY
3201	static void doSqrts(double * array, int n)
3202	{
3203	for (int i = `0`; i < n; i++)
3204	array[i] = `1.0` / sqrt(array[i]);
3205	}
3206	#endif
3207	//static int scale_stats[5]={0,0,0,0,0};
3208	// Creates scales for column copy (rowCopy in model may be modified)
3209	int
3210	ClpPackedMatrix::scale(ClpModel * model, const ClpSimplex * /baseModel/) const
3211	{
3212	//const ClpSimplex baseModel=NULL;*
3213	//return scale2(model);
3214	#if 0
3215	ClpMatrixBase * rowClone = NULL;
3216	if (model->rowCopy())
3217	rowClone = model->rowCopy()->clone();
3218	assert (!model->rowScale());
3219	assert (!model->columnScale());
3220	int returnCode = scale2(model);
3221	if (returnCode)
3222	return returnCode;
3223	#endif
3224	#ifndef NDEBUG
3225	//checkFlags();
3226	#endif
3227	int numberRows = model->numberRows();
3228	int numberColumns = matrix_->getNumCols();
3229	model->setClpScaledMatrix(NULL); // get rid of any scaled matrix
3230	// If empty - return as sanityCheck will trap
3231	if (!numberRows \|\| !numberColumns) {
3232	model->setRowScale(NULL);
3233	model->setColumnScale(NULL);
3234	return `1`;
3235	}
3236	#if 0
3237	// start fake
3238	double * rowScale2 = CoinCopyOfArray(model->rowScale(), numberRows);
3239	double * columnScale2 = CoinCopyOfArray(model->columnScale(), numberColumns);
3240	model->setRowScale(NULL);
3241	model->setColumnScale(NULL);
3242	model->setNewRowCopy(rowClone);
3243	#endif
3244	ClpMatrixBase * rowCopyBase = model->rowCopy();
3245	double * rowScale;
3246	double * columnScale;
3247	//assert (!model->rowScale());
3248	bool arraysExist;
3249	double * inverseRowScale = NULL;
3250	double * inverseColumnScale = NULL;
3251	if (!model->rowScale()) {
3252	rowScale = new double [numberRows*`2`];
3253	columnScale = new double [numberColumns*`2`];
3254	inverseRowScale = rowScale + numberRows;
3255	inverseColumnScale = columnScale + numberColumns;
3256	arraysExist = false;
3257	} else {
3258	rowScale = model->mutableRowScale();
3259	columnScale = model->mutableColumnScale();
3260	inverseRowScale = model->mutableInverseRowScale();
3261	inverseColumnScale = model->mutableInverseColumnScale();
3262	arraysExist = true;
3263	}
3264	assert (inverseRowScale == rowScale + numberRows);
3265	assert (inverseColumnScale == columnScale + numberColumns);
3266	// we are going to mark bits we are interested in
3267	char * usefulColumn = new char [numberColumns];
3268	double * rowLower = model->rowLower();
3269	double * rowUpper = model->rowUpper();
3270	double * columnLower = model->columnLower();
3271	double * columnUpper = model->columnUpper();
3272	int iColumn, iRow;
3273	//#define LEAVE_FIXED
3274	// mark empty and fixed columns
3275	// also see if worth scaling
3276	assert (model->scalingFlag() <= `5`);
3277	// scale_stats[model->scalingFlag()]++;
3278	double largest = `0.0`;
3279	double smallest = `1.0e50`;
3280	// get matrix data pointers
3281	int * row = matrix_->getMutableIndices();
3282	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3283	int * columnLength = matrix_->getMutableVectorLengths();
3284	double * elementByColumn = matrix_->getMutableElements();
3285	int deletedElements = `0`;
3286	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
3287	CoinBigIndex j;
3288	char useful = `0`;
3289	bool deleteSome = false;
3290	int start = columnStart[iColumn];
3291	int end = start + columnLength[iColumn];
3292	#ifndef LEAVE_FIXED
3293	if (columnUpper[iColumn] >
3294	columnLower[iColumn] + `1.0e-12`) {
3295	#endif
3296	for (j = start; j < end; j++) {
3297	iRow = row[j];
3298	double value = fabs(elementByColumn[j]);
3299	if (value > `1.0e-20`) {
3300	useful = `1`;
3301	largest = CoinMax(largest, fabs(elementByColumn[j]));
3302	smallest = CoinMin(smallest, fabs(elementByColumn[j]));
3303	} else {
3304	// small
3305	deleteSome = true;
3306	}
3307	}
3308	#ifndef LEAVE_FIXED
3309	} else {
3310	// just check values
3311	for (j = start; j < end; j++) {
3312	double value = fabs(elementByColumn[j]);
3313	if (value <= `1.0e-20`) {
3314	// small
3315	deleteSome = true;
3316	}
3317	}
3318	}
3319	#endif
3320	usefulColumn[iColumn] = useful;
3321	if (deleteSome) {
3322	CoinBigIndex put = start;
3323	for (j = start; j < end; j++) {
3324	double value = elementByColumn[j];
3325	if (fabs(value) > `1.0e-20`) {
3326	row[put] = row[j];
3327	elementByColumn[put++] = value;
3328	}
3329	}
3330	deletedElements += end - put;
3331	columnLength[iColumn] = put - start;
3332	}
3333	}
3334	if (deletedElements)
3335	matrix_->setNumElements(matrix_->getNumElements()-deletedElements);
3336	model->messageHandler()->message(CLP_PACKEDSCALE_INITIAL, *model->messagesPointer())
3337	<< smallest << largest
3338	<< CoinMessageEol;
3339	if (smallest >= `0.5` && largest <= `2.0` && !deletedElements) {
3340	// don't bother scaling
3341	model->messageHandler()->message(CLP_PACKEDSCALE_FORGET, *model->messagesPointer())
3342	<< CoinMessageEol;
3343	if (!arraysExist) {
3344	delete [] rowScale;
3345	delete [] columnScale;
3346	} else {
3347	model->setRowScale(NULL);
3348	model->setColumnScale(NULL);
3349	}
3350	delete [] usefulColumn;
3351	return `1`;
3352	} else {
3353	#ifdef CLP_INVESTIGATE
3354	if (deletedElements)
3355	printf("DEL_ELS\n");
3356	#endif
3357	if (!rowCopyBase) {
3358	// temporary copy
3359	rowCopyBase = reverseOrderedCopy();
3360	} else if (deletedElements) {
3361	rowCopyBase = reverseOrderedCopy();
3362	model->setNewRowCopy(rowCopyBase);
3363	}
3364	#ifndef NDEBUG
3365	ClpPackedMatrix* rowCopy =
3366	dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3367	// Make sure it is really a ClpPackedMatrix
3368	assert (rowCopy != NULL);
3369	#else
3370	ClpPackedMatrix* rowCopy =
3371	static_cast< ClpPackedMatrix*>(rowCopyBase);
3372	#endif
3373
3374	const int * column = rowCopy->getIndices();
3375	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3376	const double * element = rowCopy->getElements();
3377	// need to scale
3378	if (largest > `1.0e13` * smallest) {
3379	// safer to have smaller zero tolerance
3380	double ratio = smallest / largest;
3381	ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
3382	double newTolerance = CoinMax(ratio * `0.5`, `1.0e-18`);
3383	if (simplex->zeroTolerance() > newTolerance)
3384	simplex->setZeroTolerance(newTolerance);
3385	}
3386	int scalingMethod = model->scalingFlag();
3387	if (scalingMethod == `4`) {
3388	// As auto
3389	scalingMethod = `3`;
3390	} else if (scalingMethod == `5`) {
3391	// As geometric
3392	scalingMethod = `2`;
3393	}
3394	double savedOverallRatio = `0.0`;
3395	double tolerance = `5.0` * model->primalTolerance();
3396	double overallLargest = -`1.0e-20`;
3397	double overallSmallest = `1.0e20`;
3398	bool finished = false;
3399	// if scalingMethod 3 then may change
3400	bool extraDetails = (model->logLevel() > `2`);
3401	while (!finished) {
3402	int numberPass = `3`;
3403	overallLargest = -`1.0e-20`;
3404	overallSmallest = `1.0e20`;
3405	ClpFillN ( rowScale, numberRows, `1.0`);
3406	ClpFillN ( columnScale, numberColumns, `1.0`);
3407	if (scalingMethod == `1` \|\| scalingMethod == `3`) {
3408	// Maximum in each row
3409	for (iRow = `0`; iRow < numberRows; iRow++) {
3410	CoinBigIndex j;
3411	largest = `1.0e-10`;
3412	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
3413	int iColumn = column[j];
3414	if (usefulColumn[iColumn]) {
3415	double value = fabs(element[j]);
3416	largest = CoinMax(largest, value);
3417	assert (largest < `1.0e40`);
3418	}
3419	}
3420	rowScale[iRow] = `1.0` / largest;
3421	#ifdef COIN_DEVELOP
3422	if (extraDetails) {
3423	overallLargest = CoinMax(overallLargest, largest);
3424	overallSmallest = CoinMin(overallSmallest, largest);
3425	}
3426	#endif
3427	}
3428	} else {
3429	#ifdef USE_OBJECTIVE
3430	// This will be used to help get scale factors
3431	double * objective = new double[numberColumns];
3432	CoinMemcpyN(model->costRegion(`1`), numberColumns, objective);
3433	double objScale = `1.0`;
3434	#endif
3435	while (numberPass) {
3436	overallLargest = `0.0`;
3437	overallSmallest = `1.0e50`;
3438	numberPass--;
3439	// Geometric mean on row scales
3440	for (iRow = `0`; iRow < numberRows; iRow++) {
3441	CoinBigIndex j;
3442	largest = `1.0e-50`;
3443	smallest = `1.0e50`;
3444	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
3445	int iColumn = column[j];
3446	if (usefulColumn[iColumn]) {
3447	double value = fabs(element[j]);
3448	value *= columnScale[iColumn];
3449	largest = CoinMax(largest, value);
3450	smallest = CoinMin(smallest, value);
3451	}
3452	}
3453
3454	#ifdef SQRT_ARRAY
3455	rowScale[iRow] = smallest * largest;
3456	#else
3457	rowScale[iRow] = `1.0` / sqrt(smallest * largest);
3458	#endif
3459	//rowScale[iRow]=CoinMax(1.0e-10,CoinMin(1.0e10,rowScale[iRow]));
3460	if (extraDetails) {
3461	overallLargest = CoinMax(largest * rowScale[iRow], overallLargest);
3462	overallSmallest = CoinMin(smallest * rowScale[iRow], overallSmallest);
3463	}
3464	}
3465	if (model->scalingFlag() == `5`)
3466	break; // just scale rows
3467	#ifdef SQRT_ARRAY
3468	doSqrts(rowScale, numberRows);
3469	#endif
3470	#ifdef USE_OBJECTIVE
3471	largest = `1.0e-20`;
3472	smallest = `1.0e50`;
3473	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
3474	if (usefulColumn[iColumn]) {
3475	double value = fabs(objective[iColumn]);
3476	value *= columnScale[iColumn];
3477	largest = CoinMax(largest, value);
3478	smallest = CoinMin(smallest, value);
3479	}
3480	}
3481	objScale = `1.0` / sqrt(smallest * largest);
3482	#endif
3483	model->messageHandler()->message(CLP_PACKEDSCALE_WHILE, *model->messagesPointer())
3484	<< overallSmallest
3485	<< overallLargest
3486	<< CoinMessageEol;
3487	// skip last column round
3488	if (numberPass == `1`)
3489	break;
3490	// Geometric mean on column scales
3491	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
3492	if (usefulColumn[iColumn]) {
3493	CoinBigIndex j;
3494	largest = `1.0e-50`;
3495	smallest = `1.0e50`;
3496	for (j = columnStart[iColumn];
3497	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3498	iRow = row[j];
3499	double value = fabs(elementByColumn[j]);
3500	value *= rowScale[iRow];
3501	largest = CoinMax(largest, value);
3502	smallest = CoinMin(smallest, value);
3503	}
3504	#ifdef USE_OBJECTIVE
3505	if (fabs(objective[iColumn]) > `1.0e-20`) {
3506	double value = fabs(objective[iColumn]) * objScale;
3507	largest = CoinMax(largest, value);
3508	smallest = CoinMin(smallest, value);
3509	}
3510	#endif
3511	#ifdef SQRT_ARRAY
3512	columnScale[iColumn] = smallest * largest;
3513	#else
3514	columnScale[iColumn] = `1.0` / sqrt(smallest * largest);
3515	#endif
3516	//columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
3517	}
3518	}
3519	#ifdef SQRT_ARRAY
3520	doSqrts(columnScale, numberColumns);
3521	#endif
3522	}
3523	#ifdef USE_OBJECTIVE
3524	delete [] objective;
3525	printf("obj scale %g - use it if you want\n", objScale);
3526	#endif
3527	}
3528	// If ranges will make horrid then scale
3529	for (iRow = `0`; iRow < numberRows; iRow++) {
3530	double difference = rowUpper[iRow] - rowLower[iRow];
3531	double scaledDifference = difference * rowScale[iRow];
3532	if (scaledDifference > tolerance && scaledDifference < `1.0e-4`) {
3533	// make gap larger
3534	rowScale[iRow] *= `1.0e-4` / scaledDifference;
3535	rowScale[iRow] = CoinMax(`1.0e-10`, CoinMin(`1.0e10`, rowScale[iRow]));
3536	//printf("Row %d difference %g scaled diff %g => %g\n",iRow,difference,
3537	// scaledDifference,differencerowScale[iRow]);*
3538	}
3539	}
3540	// final pass to scale columns so largest is reasonable
3541	// See what smallest will be if largest is 1.0
3542	if (model->scalingFlag() != `5`) {
3543	overallSmallest = `1.0e50`;
3544	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
3545	if (usefulColumn[iColumn]) {
3546	CoinBigIndex j;
3547	largest = `1.0e-20`;
3548	smallest = `1.0e50`;
3549	for (j = columnStart[iColumn];
3550	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3551	iRow = row[j];
3552	double value = fabs(elementByColumn[j] * rowScale[iRow]);
3553	largest = CoinMax(largest, value);
3554	smallest = CoinMin(smallest, value);
3555	}
3556	if (overallSmallest * largest > smallest)
3557	overallSmallest = smallest / largest;
3558	}
3559	}
3560	}
3561	if (scalingMethod == `1` \|\| scalingMethod == `2`) {
3562	finished = true;
3563	} else if (savedOverallRatio == `0.0` && scalingMethod != `4`) {
3564	savedOverallRatio = overallSmallest;
3565	scalingMethod = `4`;
3566	} else {
3567	assert (scalingMethod == `4`);
3568	if (overallSmallest > `2.0` * savedOverallRatio) {
3569	finished = true; // geometric was better
3570	if (model->scalingFlag() == `4`) {
3571	// If in Branch and bound change
3572	if ((model->specialOptions() & `1024`) != `0`) {
3573	model->scaling(`2`);
3574	}
3575	}
3576	} else {
3577	scalingMethod = `1`; // redo equilibrium
3578	if (model->scalingFlag() == `4`) {
3579	// If in Branch and bound change
3580	if ((model->specialOptions() & `1024`) != `0`) {
3581	model->scaling(`1`);
3582	}
3583	}
3584	}
3585	#if 0
3586	if (extraDetails) {
3587	if (finished)
3588	printf("equilibrium ratio %g, geometric ratio %g , geo chosen\n",
3589	savedOverallRatio, overallSmallest);
3590	else
3591	printf("equilibrium ratio %g, geometric ratio %g , equi chosen\n",
3592	savedOverallRatio, overallSmallest);
3593	}
3594	#endif
3595	}
3596	}
3597	//#define RANDOMIZE
3598	#ifdef RANDOMIZE
3599	// randomize by up to 10%
3600	for (iRow = `0`; iRow < numberRows; iRow++) {
3601	double value = `0.5` - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
3602	rowScale[iRow] = (`1.0` + `0.1` value);
3603	}
3604	#endif
3605	overallLargest = `1.0`;
3606	if (overallSmallest < `1.0e-1`)
3607	overallLargest = `1.0` / sqrt(overallSmallest);
3608	overallLargest = CoinMin(`100.0`, overallLargest);
3609	overallSmallest = `1.0e50`;
3610	char * usedRow = reinterpret_cast<char *>(inverseRowScale);
3611	memset(usedRow, `0`, numberRows);
3612	//printf("scaling %d\n",model->scalingFlag());
3613	if (model->scalingFlag() != `5`) {
3614	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
3615	if (columnUpper[iColumn] >
3616	columnLower[iColumn] + `1.0e-12`) {
3617	//if (usefulColumn[iColumn]) {
3618	CoinBigIndex j;
3619	largest = `1.0e-20`;
3620	smallest = `1.0e50`;
3621	for (j = columnStart[iColumn];
3622	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
3623	iRow = row[j];
3624	usedRow[iRow] = `1`;
3625	double value = fabs(elementByColumn[j] * rowScale[iRow]);
3626	largest = CoinMax(largest, value);
3627	smallest = CoinMin(smallest, value);
3628	}
3629	columnScale[iColumn] = overallLargest / largest;
3630	//columnScale[iColumn]=CoinMax(1.0e-10,CoinMin(1.0e10,columnScale[iColumn]));
3631	#ifdef RANDOMIZE
3632	double value = `0.5` - randomNumberGenerator_.randomDouble(); //between -0.5 to + 0.5
3633	columnScale[iColumn] = (`1.0` + `0.1` value);
3634	#endif
3635	double difference = columnUpper[iColumn] - columnLower[iColumn];
3636	if (difference < `1.0e-5` * columnScale[iColumn]) {
3637	// make gap larger
3638	columnScale[iColumn] = difference / `1.0e-5`;
3639	//printf("Column %d difference %g scaled diff %g => %g\n",iColumn,difference,
3640	// scaledDifference,differencecolumnScale[iColumn]);*
3641	}
3642	double value = smallest * columnScale[iColumn];
3643	if (overallSmallest > value)
3644	overallSmallest = value;
3645	//overallSmallest = CoinMin(overallSmallest,smallestcolumnScale[iColumn]);*
3646	} else {
3647	assert(columnScale[iColumn] == `1.0`);
3648	//columnScale[iColumn]=1.0;
3649	}
3650	}
3651	for (iRow = `0`; iRow < numberRows; iRow++) {
3652	if (!usedRow[iRow]) {
3653	rowScale[iRow] = `1.0`;
3654	}
3655	}
3656	}
3657	model->messageHandler()->message(CLP_PACKEDSCALE_FINAL, *model->messagesPointer())
3658	<< overallSmallest
3659	<< overallLargest
3660	<< CoinMessageEol;
3661	#if 0
3662	{
3663	for (iRow = `0`; iRow < numberRows; iRow++) {
3664	assert (rowScale[iRow] == rowScale2[iRow]);
3665	}
3666	delete [] rowScale2;
3667	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
3668	assert (columnScale[iColumn] == columnScale2[iColumn]);
3669	}
3670	delete [] columnScale2;
3671	}
3672	#endif
3673	if (overallSmallest < `1.0e-13`) {
3674	// Change factorization zero tolerance
3675	double newTolerance = CoinMax(`1.0e-15` * (overallSmallest / `1.0e-13`),
3676	`1.0e-18`);
3677	ClpSimplex * simplex = static_cast<ClpSimplex *> (model);
3678	if (simplex->factorization()->zeroTolerance() > newTolerance)
3679	simplex->factorization()->zeroTolerance(newTolerance);
3680	newTolerance = CoinMax(overallSmallest * `0.5`, `1.0e-18`);
3681	simplex->setZeroTolerance(newTolerance);
3682	}
3683	delete [] usefulColumn;
3684	#ifndef SLIM_CLP
3685	// If quadratic then make symmetric
3686	ClpObjective * obj = model->objectiveAsObject();
3687	#ifndef NO_RTTI
3688	ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(obj));
3689	#else
3690	ClpQuadraticObjective * quadraticObj = NULL;
3691	if (obj->type() == `2`)
3692	quadraticObj = (static_cast< ClpQuadraticObjective*>(obj));
3693	#endif
3694	if (quadraticObj) {
3695	if (!rowCopyBase) {
3696	// temporary copy
3697	rowCopyBase = reverseOrderedCopy();
3698	}
3699	#ifndef NDEBUG
3700	ClpPackedMatrix* rowCopy =
3701	dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
3702	// Make sure it is really a ClpPackedMatrix
3703	assert (rowCopy != NULL);
3704	#else
3705	ClpPackedMatrix* rowCopy =
3706	static_cast< ClpPackedMatrix*>(rowCopyBase);
3707	#endif
3708	const int * column = rowCopy->getIndices();
3709	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3710	CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
3711	int numberXColumns = quadratic->getNumCols();
3712	if (numberXColumns < numberColumns) {
3713	// we assume symmetric
3714	int numberQuadraticColumns = `0`;
3715	int i;
3716	//const int columnQuadratic = quadratic->getIndices();*
3717	//const int columnQuadraticStart = quadratic->getVectorStarts();*
3718	const int * columnQuadraticLength = quadratic->getVectorLengths();
3719	for (i = `0`; i < numberXColumns; i++) {
3720	int length = columnQuadraticLength[i];
3721	#ifndef CORRECT_COLUMN_COUNTS
3722	length = `1`;
3723	#endif
3724	if (length)
3725	numberQuadraticColumns++;
3726	}
3727	int numberXRows = numberRows - numberQuadraticColumns;
3728	numberQuadraticColumns = `0`;
3729	for (i = `0`; i < numberXColumns; i++) {
3730	int length = columnQuadraticLength[i];
3731	#ifndef CORRECT_COLUMN_COUNTS
3732	length = `1`;
3733	#endif
3734	if (length) {
3735	rowScale[numberQuadraticColumns+numberXRows] = columnScale[i];
3736	numberQuadraticColumns++;
3737	}
3738	}
3739	int numberQuadraticRows = `0`;
3740	for (i = `0`; i < numberXRows; i++) {
3741	// See if any in row quadratic
3742	CoinBigIndex j;
3743	int numberQ = `0`;
3744	for (j = rowStart[i]; j < rowStart[i+`1`]; j++) {
3745	int iColumn = column[j];
3746	if (columnQuadraticLength[iColumn])
3747	numberQ++;
3748	}
3749	#ifndef CORRECT_ROW_COUNTS
3750	numberQ = `1`;
3751	#endif
3752	if (numberQ) {
3753	columnScale[numberQuadraticRows+numberXColumns] = rowScale[i];
3754	numberQuadraticRows++;
3755	}
3756	}
3757	// and make sure Sj okay
3758	for (iColumn = numberQuadraticRows + numberXColumns; iColumn < numberColumns; iColumn++) {
3759	CoinBigIndex j = columnStart[iColumn];
3760	assert(columnLength[iColumn] == `1`);
3761	int iRow = row[j];
3762	double value = fabs(elementByColumn[j] * rowScale[iRow]);
3763	columnScale[iColumn] = `1.0` / value;
3764	}
3765	}
3766	}
3767	#endif
3768	// make copy (could do faster by using previous values)
3769	// could just do partial
3770	for (iRow = `0`; iRow < numberRows; iRow++)
3771	inverseRowScale[iRow] = `1.0` / rowScale[iRow] ;
3772	for (iColumn = `0`; iColumn < numberColumns; iColumn++)
3773	inverseColumnScale[iColumn] = `1.0` / columnScale[iColumn] ;
3774	if (!arraysExist) {
3775	model->setRowScale(rowScale);
3776	model->setColumnScale(columnScale);
3777	}
3778	if (model->rowCopy()) {
3779	// need to replace row by row
3780	ClpPackedMatrix* rowCopy =
3781	static_cast< ClpPackedMatrix*>(model->rowCopy());
3782	double * element = rowCopy->getMutableElements();
3783	const int * column = rowCopy->getIndices();
3784	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
3785	// scale row copy
3786	for (iRow = `0`; iRow < numberRows; iRow++) {
3787	CoinBigIndex j;
3788	double scale = rowScale[iRow];
3789	double * elementsInThisRow = element + rowStart[iRow];
3790	const int * columnsInThisRow = column + rowStart[iRow];
3791	int number = rowStart[iRow+`1`] - rowStart[iRow];
3792	assert (number <= numberColumns);
3793	for (j = `0`; j < number; j++) {
3794	int iColumn = columnsInThisRow[j];
3795	elementsInThisRow[j] = scale columnScale[iColumn];
3796	}
3797	}
3798	if ((model->specialOptions() & `262144`) != `0`) {
3799	//if ((model->specialOptions()&(COIN_CBC_USING_CLP\|16384))!=0) {
3800	//if (model->inCbcBranchAndBound()&&false) {
3801	// copy without gaps
3802	CoinPackedMatrix * scaledMatrix = new CoinPackedMatrix (*matrix_, `0`, `0`);
3803	ClpPackedMatrix * scaled = new ClpPackedMatrix (scaledMatrix);
3804	model->setClpScaledMatrix(scaled);
3805	// get matrix data pointers
3806	const int * row = scaledMatrix->getIndices();
3807	const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
3808	#ifndef NDEBUG
3809	const int * columnLength = scaledMatrix->getVectorLengths();
3810	#endif
3811	double * elementByColumn = scaledMatrix->getMutableElements();
3812	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
3813	CoinBigIndex j;
3814	double scale = columnScale[iColumn];
3815	assert (columnStart[iColumn+`1`] == columnStart[iColumn] + columnLength[iColumn]);
3816	for (j = columnStart[iColumn];
3817	j < columnStart[iColumn+`1`]; j++) {
3818	int iRow = row[j];
3819	elementByColumn[j] = scale rowScale[iRow];
3820	}
3821	}
3822	} else {
3823	//printf("not in b&b\n");
3824	}
3825	} else {
3826	// no row copy
3827	delete rowCopyBase;
3828	}
3829	return `0`;
3830	}
3831	}
3832	// Creates scaled column copy if scales exist
3833	void
3834	ClpPackedMatrix::createScaledMatrix(ClpSimplex * model) const
3835	{
3836	int numberRows = model->numberRows();
3837	int numberColumns = matrix_->getNumCols();
3838	model->setClpScaledMatrix(NULL); // get rid of any scaled matrix
3839	// If empty - return as sanityCheck will trap
3840	if (!numberRows \|\| !numberColumns) {
3841	model->setRowScale(NULL);
3842	model->setColumnScale(NULL);
3843	return ;
3844	}
3845	if (!model->rowScale())
3846	return;
3847	double * rowScale = model->mutableRowScale();
3848	double * columnScale = model->mutableColumnScale();
3849	// copy without gaps
3850	CoinPackedMatrix * scaledMatrix = new CoinPackedMatrix (*matrix_, `0`, `0`);
3851	ClpPackedMatrix * scaled = new ClpPackedMatrix (scaledMatrix);
3852	model->setClpScaledMatrix(scaled);
3853	// get matrix data pointers
3854	const int * row = scaledMatrix->getIndices();
3855	const CoinBigIndex * columnStart = scaledMatrix->getVectorStarts();
3856	#ifndef NDEBUG
3857	const int * columnLength = scaledMatrix->getVectorLengths();
3858	#endif
3859	double * elementByColumn = scaledMatrix->getMutableElements();
3860	for (int iColumn = `0`; iColumn < numberColumns; iColumn++) {
3861	CoinBigIndex j;
3862	double scale = columnScale[iColumn];
3863	assert (columnStart[iColumn+`1`] == columnStart[iColumn] + columnLength[iColumn]);
3864	for (j = columnStart[iColumn];
3865	j < columnStart[iColumn+`1`]; j++) {
3866	int iRow = row[j];
3867	elementByColumn[j] = scale rowScale[iRow];
3868	}
3869	}
3870	}
3871	/ Unpacks a column into an CoinIndexedvector*
3872	*/
3873	void
3874	ClpPackedMatrix::unpack(const ClpSimplex * model, CoinIndexedVector * rowArray,
3875	int iColumn) const
3876	{
3877	const double * rowScale = model->rowScale();
3878	const int * row = matrix_->getIndices();
3879	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3880	const int * columnLength = matrix_->getVectorLengths();
3881	const double * elementByColumn = matrix_->getElements();
3882	CoinBigIndex i;
3883	if (!rowScale) {
3884	for (i = columnStart[iColumn];
3885	i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3886	rowArray->add(row[i], elementByColumn[i]);
3887	}
3888	} else {
3889	// apply scaling
3890	double scale = model->columnScale()[iColumn];
3891	for (i = columnStart[iColumn];
3892	i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3893	int iRow = row[i];
3894	rowArray->add(iRow, elementByColumn[i]scale rowScale[iRow]);
3895	}
3896	}
3897	}
3898	/ Unpacks a column into a CoinIndexedVector*
3899	** in packed format
3900	Note that model is NOT const. Bounds and objective could
3901	be modified if doing column generation (just for this variable) /*
3902	void
3903	ClpPackedMatrix::unpackPacked(ClpSimplex * model,
3904	CoinIndexedVector * rowArray,
3905	int iColumn) const
3906	{
3907	const double * rowScale = model->rowScale();
3908	const int * row = matrix_->getIndices();
3909	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3910	const int * columnLength = matrix_->getVectorLengths();
3911	const double * elementByColumn = matrix_->getElements();
3912	CoinBigIndex i;
3913	int * index = rowArray->getIndices();
3914	double * array = rowArray->denseVector();
3915	int number = `0`;
3916	if (!rowScale) {
3917	for (i = columnStart[iColumn];
3918	i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3919	int iRow = row[i];
3920	double value = elementByColumn[i];
3921	if (value) {
3922	array[number] = value;
3923	index[number++] = iRow;
3924	}
3925	}
3926	rowArray->setNumElements(number);
3927	rowArray->setPackedMode(true);
3928	} else {
3929	// apply scaling
3930	double scale = model->columnScale()[iColumn];
3931	for (i = columnStart[iColumn];
3932	i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3933	int iRow = row[i];
3934	double value = elementByColumn[i] * scale * rowScale[iRow];
3935	if (value) {
3936	array[number] = value;
3937	index[number++] = iRow;
3938	}
3939	}
3940	rowArray->setNumElements(number);
3941	rowArray->setPackedMode(true);
3942	}
3943	}
3944	/ Adds multiple of a column into an CoinIndexedvector*
3945	You can use quickAdd to add to vector /*
3946	void
3947	ClpPackedMatrix::add(const ClpSimplex * model, CoinIndexedVector * rowArray,
3948	int iColumn, double multiplier) const
3949	{
3950	const double * rowScale = model->rowScale();
3951	const int * row = matrix_->getIndices();
3952	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3953	const int * columnLength = matrix_->getVectorLengths();
3954	const double * elementByColumn = matrix_->getElements();
3955	CoinBigIndex i;
3956	if (!rowScale) {
3957	for (i = columnStart[iColumn];
3958	i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3959	int iRow = row[i];
3960	rowArray->quickAdd(iRow, multiplier * elementByColumn[i]);
3961	}
3962	} else {
3963	// apply scaling
3964	double scale = model->columnScale()[iColumn] * multiplier;
3965	for (i = columnStart[iColumn];
3966	i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3967	int iRow = row[i];
3968	rowArray->quickAdd(iRow, elementByColumn[i]scale rowScale[iRow]);
3969	}
3970	}
3971	}
3972	/ Adds multiple of a column into an array /
3973	void
3974	ClpPackedMatrix::add(const ClpSimplex * model, double * array,
3975	int iColumn, double multiplier) const
3976	{
3977	const double * rowScale = model->rowScale();
3978	const int * row = matrix_->getIndices();
3979	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
3980	const int * columnLength = matrix_->getVectorLengths();
3981	const double * elementByColumn = matrix_->getElements();
3982	CoinBigIndex i;
3983	if (!rowScale) {
3984	for (i = columnStart[iColumn];
3985	i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3986	int iRow = row[i];
3987	array[iRow] += multiplier * elementByColumn[i];
3988	}
3989	} else {
3990	// apply scaling
3991	double scale = model->columnScale()[iColumn] * multiplier;
3992	for (i = columnStart[iColumn];
3993	i < columnStart[iColumn] + columnLength[iColumn]; i++) {
3994	int iRow = row[i];
3995	array[iRow] += elementByColumn[i] * scale * rowScale[iRow];
3996	}
3997	}
3998	}
3999	/ Checks if all elements are in valid range. Can just*
4000	return true if you are not paranoid. For Clp I will
4001	probably expect no zeros. Code can modify matrix to get rid of
4002	small elements.
4003	check bits (can be turned off to save time) :
4004	1 - check if matrix has gaps
4005	2 - check if zero elements
4006	4 - check and compress duplicates
4007	8 - report on large and small
4008	*/
4009	bool
4010	ClpPackedMatrix::allElementsInRange(ClpModel * model,
4011	double smallest, double largest,
4012	int check)
4013	{
4014	int iColumn;
4015	// make sure matrix correct size
4016	assert (matrix_->getNumRows() <= model->numberRows());
4017	if (model->clpScaledMatrix())
4018	assert (model->clpScaledMatrix()->getNumElements() == matrix_->getNumElements());
4019	assert (matrix_->getNumRows() <= model->numberRows());
4020	matrix_->setDimensions(model->numberRows(), model->numberColumns());
4021	CoinBigIndex numberLarge = `0`;
4022	CoinBigIndex numberSmall = `0`;
4023	CoinBigIndex numberDuplicate = `0`;
4024	int firstBadColumn = -`1`;
4025	int firstBadRow = -`1`;
4026	double firstBadElement = `0.0`;
4027	// get matrix data pointers
4028	const int * row = matrix_->getIndices();
4029	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4030	const int * columnLength = matrix_->getVectorLengths();
4031	const double * elementByColumn = matrix_->getElements();
4032	int numberRows = model->numberRows();
4033	int numberColumns = matrix_->getNumCols();
4034	// Say no gaps
4035	flags_ &= ~`2`;
4036	if (type_>=`10`)
4037	return true; // gub
4038	if (check == `14` \|\| check == `10`) {
4039	if (matrix_->getNumElements() < columnStart[numberColumns]) {
4040	// pack down!
4041	#if 0
4042	matrix_->removeGaps();
4043	#else
4044	checkGaps();
4045	#endif
4046	#ifdef COIN_DEVELOP
4047	//printf("flags set to 2\n");
4048	#endif
4049	} else if (numberColumns) {
4050	assert(columnStart[numberColumns] ==
4051	columnStart[numberColumns-`1`] + columnLength[numberColumns-`1`]);
4052	}
4053	return true;
4054	}
4055	assert (check == `15` \|\| check == `11`);
4056	if (check == `15`) {
4057	int * mark = new int [numberRows];
4058	int i;
4059	for (i = `0`; i < numberRows; i++)
4060	mark[i] = -`1`;
4061	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
4062	CoinBigIndex j;
4063	CoinBigIndex start = columnStart[iColumn];
4064	CoinBigIndex end = start + columnLength[iColumn];
4065	if (end != columnStart[iColumn+`1`])
4066	flags_ \|= `2`;
4067	for (j = start; j < end; j++) {
4068	double value = fabs(elementByColumn[j]);
4069	int iRow = row[j];
4070	if (iRow < `0` \|\| iRow >= numberRows) {
4071	#ifndef COIN_BIG_INDEX
4072	printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4073	#elif COIN_BIG_INDEX==0
4074	printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4075	#elif COIN_BIG_INDEX==1
4076	printf("Out of range %d %ld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4077	#else
4078	printf("Out of range %d %lld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4079	#endif
4080	delete [] mark;
4081	return false;
4082	}
4083	if (mark[iRow] == -`1`) {
4084	mark[iRow] = j;
4085	} else {
4086	// duplicate
4087	numberDuplicate++;
4088	}
4089	//printf("%d %d %d %g\n",iColumn,j,row[j],elementByColumn[j]);
4090	if (!value)
4091	flags_ \|= `1`; // there are zero elements
4092	if (value < smallest) {
4093	numberSmall++;
4094	} else if (!(value <= largest)) {
4095	numberLarge++;
4096	if (firstBadColumn < `0`) {
4097	firstBadColumn = iColumn;
4098	firstBadRow = row[j];
4099	firstBadElement = elementByColumn[j];
4100	}
4101	}
4102	}
4103	//clear mark
4104	for (j = columnStart[iColumn];
4105	j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4106	int iRow = row[j];
4107	mark[iRow] = -`1`;
4108	}
4109	}
4110	delete [] mark;
4111	} else {
4112	// just check for out of range - not for duplicates
4113	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
4114	CoinBigIndex j;
4115	CoinBigIndex start = columnStart[iColumn];
4116	CoinBigIndex end = start + columnLength[iColumn];
4117	if (end != columnStart[iColumn+`1`])
4118	flags_ \|= `2`;
4119	for (j = start; j < end; j++) {
4120	double value = fabs(elementByColumn[j]);
4121	int iRow = row[j];
4122	if (iRow < `0` \|\| iRow >= numberRows) {
4123	#ifndef COIN_BIG_INDEX
4124	printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4125	#elif COIN_BIG_INDEX==0
4126	printf("Out of range %d %d %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4127	#elif COIN_BIG_INDEX==1
4128	printf("Out of range %d %ld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4129	#else
4130	printf("Out of range %d %lld %d %g\n", iColumn, j, row[j], elementByColumn[j]);
4131	#endif
4132	return false;
4133	}
4134	if (!value)
4135	flags_ \|= `1`; // there are zero elements
4136	if (value < smallest) {
4137	numberSmall++;
4138	} else if (!(value <= largest)) {
4139	numberLarge++;
4140	if (firstBadColumn < `0`) {
4141	firstBadColumn = iColumn;
4142	firstBadRow = iRow;
4143	firstBadElement = value;
4144	}
4145	}
4146	}
4147	}
4148	}
4149	if (numberLarge) {
4150	model->messageHandler()->message(CLP_BAD_MATRIX, model->messages())
4151	<< numberLarge
4152	<< firstBadColumn << firstBadRow << firstBadElement
4153	<< CoinMessageEol;
4154	return false;
4155	}
4156	if (numberSmall)
4157	model->messageHandler()->message(CLP_SMALLELEMENTS, model->messages())
4158	<< numberSmall
4159	<< CoinMessageEol;
4160	if (numberDuplicate)
4161	model->messageHandler()->message(CLP_DUPLICATEELEMENTS, model->messages())
4162	<< numberDuplicate
4163	<< CoinMessageEol;
4164	if (numberDuplicate)
4165	matrix_->eliminateDuplicates(smallest);
4166	else if (numberSmall)
4167	matrix_->compress(smallest);
4168	// If smallest >0.0 then there can't be zero elements
4169	if (smallest > `0.0`)
4170	flags_ &= ~`1`;
4171	if (numberSmall \|\| numberDuplicate)
4172	flags_ \|= `2`; // will have gaps
4173	return true;
4174	}
4175	int
4176	ClpPackedMatrix::gutsOfTransposeTimesByRowGE3a(const CoinIndexedVector * COIN_RESTRICT piVector,
4177	int * COIN_RESTRICT index,
4178	double * COIN_RESTRICT output,
4179	int * COIN_RESTRICT lookup,
4180	char * COIN_RESTRICT marked,
4181	const double tolerance,
4182	const double scalar) const
4183	{
4184	const double * COIN_RESTRICT pi = piVector->denseVector();
4185	int numberNonZero = `0`;
4186	int numberInRowArray = piVector->getNumElements();
4187	const int * COIN_RESTRICT column = matrix_->getIndices();
4188	const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
4189	const double * COIN_RESTRICT element = matrix_->getElements();
4190	const int * COIN_RESTRICT whichRow = piVector->getIndices();
4191	int * fakeRow = const_cast<int *> (whichRow);
4192	fakeRow[numberInRowArray]=`0`; // so can touch
4193	#ifndef NDEBUG
4194	int maxColumn = `0`;
4195	#endif
4196	// * Row copy is already scaled*
4197	int nextRow=whichRow[`0`];
4198	CoinBigIndex nextStart = rowStart[nextRow];
4199	CoinBigIndex nextEnd = rowStart[nextRow+`1`];
4200	for (int i = `0`; i < numberInRowArray; i++) {
4201	double value = pi[i] * scalar;
4202	CoinBigIndex start=nextStart;
4203	CoinBigIndex end=nextEnd;
4204	nextRow=whichRow[i+`1`];
4205	nextStart = rowStart[nextRow];
4206	//coin_prefetch_const(column + nextStart);
4207	//coin_prefetch_const(element + nextStart);
4208	nextEnd = rowStart[nextRow+`1`];
4209	CoinBigIndex j;
4210	for (j = start; j < end; j++) {
4211	int iColumn = column[j];
4212	#ifndef NDEBUG
4213	maxColumn = CoinMax(maxColumn, iColumn);
4214	#endif
4215	double elValue = element[j];
4216	elValue *= value;
4217	if (marked[iColumn]) {
4218	int k = lookup[iColumn];
4219	output[k] += elValue;
4220	} else {
4221	output[numberNonZero] = elValue;
4222	marked[iColumn] = `1`;
4223	lookup[iColumn] = numberNonZero;
4224	index[numberNonZero++] = iColumn;
4225	}
4226	}
4227	}
4228	#ifndef NDEBUG
4229	int saveN = numberNonZero;
4230	#endif
4231	// get rid of tiny values and zero out lookup
4232	for (int i = `0`; i < numberNonZero; i++) {
4233	int iColumn = index[i];
4234	marked[iColumn] = `0`;
4235	double value = output[i];
4236	if (fabs(value) <= tolerance) {
4237	while (fabs(value) <= tolerance) {
4238	numberNonZero--;
4239	value = output[numberNonZero];
4240	iColumn = index[numberNonZero];
4241	marked[iColumn] = `0`;
4242	if (i < numberNonZero) {
4243	output[numberNonZero] = `0.0`;
4244	output[i] = value;
4245	index[i] = iColumn;
4246	} else {
4247	output[i] = `0.0`;
4248	value = `1.0`; // to force end of while
4249	}
4250	}
4251	}
4252	}
4253	#ifndef NDEBUG
4254	for (int i = numberNonZero; i < saveN; i++)
4255	assert(!output[i]);
4256	for (int i = `0`; i <= maxColumn; i++)
4257	assert (!marked[i]);
4258	#endif
4259	return numberNonZero;
4260	}
4261	int
4262	ClpPackedMatrix::gutsOfTransposeTimesByRowGE3(const CoinIndexedVector * COIN_RESTRICT piVector,
4263	int * COIN_RESTRICT index,
4264	double * COIN_RESTRICT output,
4265	double * COIN_RESTRICT array,
4266	const double tolerance,
4267	const double scalar) const
4268	{
4269	const double * COIN_RESTRICT pi = piVector->denseVector();
4270	int numberNonZero = `0`;
4271	int numberInRowArray = piVector->getNumElements();
4272	const int * COIN_RESTRICT column = matrix_->getIndices();
4273	const CoinBigIndex * COIN_RESTRICT rowStart = matrix_->getVectorStarts();
4274	const double * COIN_RESTRICT element = matrix_->getElements();
4275	const int * COIN_RESTRICT whichRow = piVector->getIndices();
4276	// * Row copy is already scaled*
4277	for (int i = `0`; i < numberInRowArray; i++) {
4278	int iRow = whichRow[i];
4279	double value = pi[i] * scalar;
4280	CoinBigIndex j;
4281	for (j = rowStart[iRow]; j < rowStart[iRow+`1`]; j++) {
4282	int iColumn = column[j];
4283	double inValue = array[iColumn];
4284	double elValue = element[j];
4285	elValue *= value;
4286	if (inValue) {
4287	double outValue = inValue + elValue;
4288	if (!outValue)
4289	outValue = COIN_INDEXED_REALLY_TINY_ELEMENT;
4290	array[iColumn] = outValue;
4291	} else {
4292	array[iColumn] = elValue;
4293	assert (elValue);
4294	index[numberNonZero++] = iColumn;
4295	}
4296	}
4297	}
4298	int saveN = numberNonZero;
4299	// get rid of tiny values
4300	numberNonZero=`0`;
4301	for (int i = `0`; i < saveN; i++) {
4302	int iColumn = index[i];
4303	double value = array[iColumn];
4304	array[iColumn] =`0.0`;
4305	if (fabs(value) > tolerance) {
4306	output[numberNonZero] = value;
4307	index[numberNonZero++] = iColumn;
4308	}
4309	}
4310	return numberNonZero;
4311	}
4312	/ Given positive integer weights for each row fills in sum of weights*
4313	for each column (and slack).
4314	Returns weights vector
4315	*/
4316	CoinBigIndex *
4317	ClpPackedMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
4318	{
4319	int numberRows = model->numberRows();
4320	int numberColumns = matrix_->getNumCols();
4321	int number = numberRows + numberColumns;
4322	CoinBigIndex * weights = new CoinBigIndex[number];
4323	// get matrix data pointers
4324	const int * row = matrix_->getIndices();
4325	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4326	const int * columnLength = matrix_->getVectorLengths();
4327	int i;
4328	for (i = `0`; i < numberColumns; i++) {
4329	CoinBigIndex j;
4330	CoinBigIndex count = `0`;
4331	for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
4332	int iRow = row[j];
4333	count += inputWeights[iRow];
4334	}
4335	weights[i] = count;
4336	}
4337	for (i = `0`; i < numberRows; i++) {
4338	weights[i+numberColumns] = inputWeights[i];
4339	}
4340	return weights;
4341	}
4342	/ Returns largest and smallest elements of both signs.*
4343	Largest refers to largest absolute value.
4344	*/
4345	void
4346	ClpPackedMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
4347	double & smallestPositive, double & largestPositive)
4348	{
4349	smallestNegative = -COIN_DBL_MAX;
4350	largestNegative = `0.0`;
4351	smallestPositive = COIN_DBL_MAX;
4352	largestPositive = `0.0`;
4353	// get matrix data pointers
4354	const double * elementByColumn = matrix_->getElements();
4355	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4356	const int * columnLength = matrix_->getVectorLengths();
4357	int numberColumns = matrix_->getNumCols();
4358	int i;
4359	for (i = `0`; i < numberColumns; i++) {
4360	CoinBigIndex j;
4361	for (j = columnStart[i]; j < columnStart[i] + columnLength[i]; j++) {
4362	double value = elementByColumn[j];
4363	if (value > `0.0`) {
4364	smallestPositive = CoinMin(smallestPositive, value);
4365	largestPositive = CoinMax(largestPositive, value);
4366	} else if (value < `0.0`) {
4367	smallestNegative = CoinMax(smallestNegative, value);
4368	largestNegative = CoinMin(largestNegative, value);
4369	}
4370	}
4371	}
4372	}
4373	// Says whether it can do partial pricing
4374	bool
4375	ClpPackedMatrix::canDoPartialPricing() const
4376	{
4377	return true;
4378	}
4379	// Partial pricing
4380	void
4381	ClpPackedMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
4382	int & bestSequence, int & numberWanted)
4383	{
4384	numberWanted = currentWanted_;
4385	int start = static_cast<int> (startFraction * numberActiveColumns_);
4386	int end = CoinMin(static_cast<int> (endFraction * numberActiveColumns_ + `1`), numberActiveColumns_);
4387	const double * element = matrix_->getElements();
4388	const int * row = matrix_->getIndices();
4389	const CoinBigIndex * startColumn = matrix_->getVectorStarts();
4390	const int * length = matrix_->getVectorLengths();
4391	const double * rowScale = model->rowScale();
4392	const double * columnScale = model->columnScale();
4393	int iSequence;
4394	CoinBigIndex j;
4395	double tolerance = model->currentDualTolerance();
4396	double * reducedCost = model->djRegion();
4397	const double * duals = model->dualRowSolution();
4398	const double * cost = model->costRegion();
4399	double bestDj;
4400	if (bestSequence >= `0`)
4401	bestDj = fabs(model->clpMatrix()->reducedCost(model, bestSequence));
4402	else
4403	bestDj = tolerance;
4404	int sequenceOut = model->sequenceOut();
4405	int saveSequence = bestSequence;
4406	int lastScan = minimumObjectsScan_ < `0` ? end : start + minimumObjectsScan_;
4407	int minNeg = minimumGoodReducedCosts_ == -`1` ? numberWanted : minimumGoodReducedCosts_;
4408	if (rowScale) {
4409	// scaled
4410	for (iSequence = start; iSequence < end; iSequence++) {
4411	if (iSequence != sequenceOut) {
4412	double value;
4413	ClpSimplex::Status status = model->getStatus(iSequence);
4414
4415	switch(status) {
4416
4417	case ClpSimplex::basic:
4418	case ClpSimplex::isFixed:
4419	break;
4420	case ClpSimplex::isFree:
4421	case ClpSimplex::superBasic:
4422	value = `0.0`;
4423	// scaled
4424	for (j = startColumn[iSequence];
4425	j < startColumn[iSequence] + length[iSequence]; j++) {
4426	int jRow = row[j];
4427	value -= duals[jRow] * element[j] * rowScale[jRow];
4428	}
4429	value = fabs(cost[iSequence] + value * columnScale[iSequence]);
4430	if (value > FREE_ACCEPT * tolerance) {
4431	numberWanted--;
4432	// we are going to bias towards free (but only if reasonable)
4433	value *= FREE_BIAS;
4434	if (value > bestDj) {
4435	// check flagged variable and correct dj
4436	if (!model->flagged(iSequence)) {
4437	bestDj = value;
4438	bestSequence = iSequence;
4439	} else {
4440	// just to make sure we don't exit before got something
4441	numberWanted++;
4442	}
4443	}
4444	}
4445	break;
4446	case ClpSimplex::atUpperBound:
4447	value = `0.0`;
4448	// scaled
4449	for (j = startColumn[iSequence];
4450	j < startColumn[iSequence] + length[iSequence]; j++) {
4451	int jRow = row[j];
4452	value -= duals[jRow] * element[j] * rowScale[jRow];
4453	}
4454	value = cost[iSequence] + value * columnScale[iSequence];
4455	if (value > tolerance) {
4456	numberWanted--;
4457	if (value > bestDj) {
4458	// check flagged variable and correct dj
4459	if (!model->flagged(iSequence)) {
4460	bestDj = value;
4461	bestSequence = iSequence;
4462	} else {
4463	// just to make sure we don't exit before got something
4464	numberWanted++;
4465	}
4466	}
4467	}
4468	break;
4469	case ClpSimplex::atLowerBound:
4470	value = `0.0`;
4471	// scaled
4472	for (j = startColumn[iSequence];
4473	j < startColumn[iSequence] + length[iSequence]; j++) {
4474	int jRow = row[j];
4475	value -= duals[jRow] * element[j] * rowScale[jRow];
4476	}
4477	value = -(cost[iSequence] + value * columnScale[iSequence]);
4478	if (value > tolerance) {
4479	numberWanted--;
4480	if (value > bestDj) {
4481	// check flagged variable and correct dj
4482	if (!model->flagged(iSequence)) {
4483	bestDj = value;
4484	bestSequence = iSequence;
4485	} else {
4486	// just to make sure we don't exit before got something
4487	numberWanted++;
4488	}
4489	}
4490	}
4491	break;
4492	}
4493	}
4494	if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
4495	// give up
4496	break;
4497	}
4498	if (!numberWanted)
4499	break;
4500	}
4501	if (bestSequence != saveSequence) {
4502	// recompute dj
4503	double value = `0.0`;
4504	// scaled
4505	for (j = startColumn[bestSequence];
4506	j < startColumn[bestSequence] + length[bestSequence]; j++) {
4507	int jRow = row[j];
4508	value -= duals[jRow] * element[j] * rowScale[jRow];
4509	}
4510	reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence];
4511	savedBestSequence_ = bestSequence;
4512	savedBestDj_ = reducedCost[savedBestSequence_];
4513	}
4514	} else {
4515	// not scaled
4516	for (iSequence = start; iSequence < end; iSequence++) {
4517	if (iSequence != sequenceOut) {
4518	double value;
4519	ClpSimplex::Status status = model->getStatus(iSequence);
4520
4521	switch(status) {
4522
4523	case ClpSimplex::basic:
4524	case ClpSimplex::isFixed:
4525	break;
4526	case ClpSimplex::isFree:
4527	case ClpSimplex::superBasic:
4528	value = cost[iSequence];
4529	for (j = startColumn[iSequence];
4530	j < startColumn[iSequence] + length[iSequence]; j++) {
4531	int jRow = row[j];
4532	value -= duals[jRow] * element[j];
4533	}
4534	value = fabs(value);
4535	if (value > FREE_ACCEPT * tolerance) {
4536	numberWanted--;
4537	// we are going to bias towards free (but only if reasonable)
4538	value *= FREE_BIAS;
4539	if (value > bestDj) {
4540	// check flagged variable and correct dj
4541	if (!model->flagged(iSequence)) {
4542	bestDj = value;
4543	bestSequence = iSequence;
4544	} else {
4545	// just to make sure we don't exit before got something
4546	numberWanted++;
4547	}
4548	}
4549	}
4550	break;
4551	case ClpSimplex::atUpperBound:
4552	value = cost[iSequence];
4553	// scaled
4554	for (j = startColumn[iSequence];
4555	j < startColumn[iSequence] + length[iSequence]; j++) {
4556	int jRow = row[j];
4557	value -= duals[jRow] * element[j];
4558	}
4559	if (value > tolerance) {
4560	numberWanted--;
4561	if (value > bestDj) {
4562	// check flagged variable and correct dj
4563	if (!model->flagged(iSequence)) {
4564	bestDj = value;
4565	bestSequence = iSequence;
4566	} else {
4567	// just to make sure we don't exit before got something
4568	numberWanted++;
4569	}
4570	}
4571	}
4572	break;
4573	case ClpSimplex::atLowerBound:
4574	value = cost[iSequence];
4575	for (j = startColumn[iSequence];
4576	j < startColumn[iSequence] + length[iSequence]; j++) {
4577	int jRow = row[j];
4578	value -= duals[jRow] * element[j];
4579	}
4580	value = -value;
4581	if (value > tolerance) {
4582	numberWanted--;
4583	if (value > bestDj) {
4584	// check flagged variable and correct dj
4585	if (!model->flagged(iSequence)) {
4586	bestDj = value;
4587	bestSequence = iSequence;
4588	} else {
4589	// just to make sure we don't exit before got something
4590	numberWanted++;
4591	}
4592	}
4593	}
4594	break;
4595	}
4596	}
4597	if (numberWanted + minNeg < originalWanted_ && iSequence > lastScan) {
4598	// give up
4599	break;
4600	}
4601	if (!numberWanted)
4602	break;
4603	}
4604	if (bestSequence != saveSequence) {
4605	// recompute dj
4606	double value = cost[bestSequence];
4607	for (j = startColumn[bestSequence];
4608	j < startColumn[bestSequence] + length[bestSequence]; j++) {
4609	int jRow = row[j];
4610	value -= duals[jRow] * element[j];
4611	}
4612	reducedCost[bestSequence] = value;
4613	savedBestSequence_ = bestSequence;
4614	savedBestDj_ = reducedCost[savedBestSequence_];
4615	}
4616	}
4617	currentWanted_ = numberWanted;
4618	}
4619	// Sets up an effective RHS
4620	void
4621	ClpPackedMatrix::useEffectiveRhs(ClpSimplex * model)
4622	{
4623	delete [] rhsOffset_;
4624	int numberRows = model->numberRows();
4625	rhsOffset_ = new double[numberRows];
4626	rhsOffset(model, true);
4627	}
4628	// Gets rid of special copies
4629	void
4630	ClpPackedMatrix::clearCopies()
4631	{
4632	delete rowCopy_;
4633	delete columnCopy_;
4634	rowCopy_ = NULL;
4635	columnCopy_ = NULL;
4636	flags_ &= ~(`4` + `8`);
4637	checkGaps();
4638	}
4639	// makes sure active columns correct
4640	int
4641	ClpPackedMatrix::refresh(ClpSimplex * )
4642	{
4643	numberActiveColumns_ = matrix_->getNumCols();
4644	#if 0
4645	ClpMatrixBase * rowCopyBase = reverseOrderedCopy();
4646	ClpPackedMatrix* rowCopy =
4647	dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
4648	// Make sure it is really a ClpPackedMatrix
4649	assert (rowCopy != NULL);
4650
4651	const int * column = rowCopy->matrix_->getIndices();
4652	const CoinBigIndex * rowStart = rowCopy->matrix_->getVectorStarts();
4653	const int * rowLength = rowCopy->matrix_->getVectorLengths();
4654	const double * element = rowCopy->matrix_->getElements();
4655	int numberRows = rowCopy->matrix_->getNumRows();
4656	for (int i = `0`; i < numberRows; i++) {
4657	if (!rowLength[i])
4658	printf("zero row %d\n", i);
4659	}
4660	delete rowCopy;
4661	#endif
4662	return `0`;
4663	}
4664
4665	/ Scales rowCopy if column copy scaled*
4666	Only called if scales already exist /*
4667	void
4668	ClpPackedMatrix::scaleRowCopy(ClpModel * model) const
4669	{
4670	if (model->rowCopy()) {
4671	// need to replace row by row
4672	int numberRows = model->numberRows();
4673	#ifndef NDEBUG
4674	int numberColumns = matrix_->getNumCols();
4675	#endif
4676	ClpMatrixBase * rowCopyBase = model->rowCopy();
4677	#ifndef NDEBUG
4678	ClpPackedMatrix* rowCopy =
4679	dynamic_cast< ClpPackedMatrix*>(rowCopyBase);
4680	// Make sure it is really a ClpPackedMatrix
4681	assert (rowCopy != NULL);
4682	#else
4683	ClpPackedMatrix* rowCopy =
4684	static_cast< ClpPackedMatrix*>(rowCopyBase);
4685	#endif
4686
4687	const int * column = rowCopy->getIndices();
4688	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
4689	double * element = rowCopy->getMutableElements();
4690	const double * rowScale = model->rowScale();
4691	const double * columnScale = model->columnScale();
4692	// scale row copy
4693	for (int iRow = `0`; iRow < numberRows; iRow++) {
4694	CoinBigIndex j;
4695	double scale = rowScale[iRow];
4696	double * elementsInThisRow = element + rowStart[iRow];
4697	const int * columnsInThisRow = column + rowStart[iRow];
4698	int number = rowStart[iRow+`1`] - rowStart[iRow];
4699	assert (number <= numberColumns);
4700	for (j = `0`; j < number; j++) {
4701	int iColumn = columnsInThisRow[j];
4702	elementsInThisRow[j] = scale columnScale[iColumn];
4703	}
4704	}
4705	}
4706	}
4707	/ Realy really scales column copy*
4708	Only called if scales already exist.
4709	Up to user ro delete /*
4710	ClpMatrixBase *
4711	ClpPackedMatrix::scaledColumnCopy(ClpModel * model) const
4712	{
4713	// need to replace column by column
4714	#ifndef NDEBUG
4715	int numberRows = model->numberRows();
4716	#endif
4717	int numberColumns = matrix_->getNumCols();
4718	ClpPackedMatrix * copy = new ClpPackedMatrix (*this);
4719	const int * row = copy->getIndices();
4720	const CoinBigIndex * columnStart = copy->getVectorStarts();
4721	const int * length = copy->getVectorLengths();
4722	double * element = copy->getMutableElements();
4723	const double * rowScale = model->rowScale();
4724	const double * columnScale = model->columnScale();
4725	// scale column copy
4726	for (int iColumn = `0`; iColumn < numberColumns; iColumn++) {
4727	CoinBigIndex j;
4728	double scale = columnScale[iColumn];
4729	double * elementsInThisColumn = element + columnStart[iColumn];
4730	const int * rowsInThisColumn = row + columnStart[iColumn];
4731	int number = length[iColumn];
4732	assert (number <= numberRows);
4733	for (j = `0`; j < number; j++) {
4734	int iRow = rowsInThisColumn[j];
4735	elementsInThisColumn[j] = scale rowScale[iRow];
4736	}
4737	}
4738	return copy;
4739	}
4740	// Really scale matrix
4741	void
4742	ClpPackedMatrix::reallyScale(const double * rowScale, const double * columnScale)
4743	{
4744	clearCopies();
4745	int numberColumns = matrix_->getNumCols();
4746	const int * row = matrix_->getIndices();
4747	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4748	const int * length = matrix_->getVectorLengths();
4749	double * element = matrix_->getMutableElements();
4750	// scale
4751	for (int iColumn = `0`; iColumn < numberColumns; iColumn++) {
4752	CoinBigIndex j;
4753	double scale = columnScale[iColumn];
4754	for (j = columnStart[iColumn]; j < columnStart[iColumn] + length[iColumn]; j++) {
4755	int iRow = row[j];
4756	element[j] = scale rowScale[iRow];
4757	}
4758	}
4759	}
4760	/ Delete the columns whose indices are listed in <code>indDel</code>. /
4761	void
4762	ClpPackedMatrix::deleteCols(const int numDel, const int * indDel)
4763	{
4764	if (matrix_->getNumCols())
4765	matrix_->deleteCols(numDel, indDel);
4766	clearCopies();
4767	numberActiveColumns_ = matrix_->getNumCols();
4768	// may now have gaps
4769	checkGaps();
4770	matrix_->setExtraGap(`0.0`);
4771	}
4772	/ Delete the rows whose indices are listed in <code>indDel</code>. /
4773	void
4774	ClpPackedMatrix::deleteRows(const int numDel, const int * indDel)
4775	{
4776	if (matrix_->getNumRows())
4777	matrix_->deleteRows(numDel, indDel);
4778	clearCopies();
4779	numberActiveColumns_ = matrix_->getNumCols();
4780	// may now have gaps
4781	checkGaps();
4782	matrix_->setExtraGap(`0.0`);
4783	}
4784	#ifndef CLP_NO_VECTOR
4785	// Append Columns
4786	void
4787	ClpPackedMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
4788	{
4789	matrix_->appendCols(number, columns);
4790	numberActiveColumns_ = matrix_->getNumCols();
4791	clearCopies();
4792	}
4793	// Append Rows
4794	void
4795	ClpPackedMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
4796	{
4797	matrix_->appendRows(number, rows);
4798	numberActiveColumns_ = matrix_->getNumCols();
4799	// may now have gaps
4800	checkGaps();
4801	clearCopies();
4802	}
4803	#endif
4804	/ Set the dimensions of the matrix. In effect, append new empty*
4805	columns/rows to the matrix. A negative number for either dimension
4806	means that that dimension doesn't change. Otherwise the new dimensions
4807	MUST be at least as large as the current ones otherwise an exception
4808	is thrown. /*
4809	void
4810	ClpPackedMatrix::setDimensions(int numrows, int numcols)
4811	{
4812	matrix_->setDimensions(numrows, numcols);
4813	}
4814	/ Append a set of rows/columns to the end of the matrix. Returns number of errors*
4815	i.e. if any of the new rows/columns contain an index that's larger than the
4816	number of columns-1/rows-1 (if numberOther>0) or duplicates
4817	If 0 then rows, 1 if columns /*
4818	int
4819	ClpPackedMatrix::appendMatrix(int number, int type,
4820	const CoinBigIndex * starts, const int * index,
4821	const double * element, int numberOther)
4822	{
4823	int numberErrors = `0`;
4824	// make sure other dimension is big enough
4825	if (type == `0`) {
4826	// rows
4827	if (matrix_->isColOrdered() && numberOther > matrix_->getNumCols())
4828	matrix_->setDimensions(-`1`, numberOther);
4829	if (!matrix_->isColOrdered() \|\| numberOther >= `0` \|\| matrix_->getExtraGap() \|\| matrix_->hasGaps()) {
4830	numberErrors = matrix_->appendRows(number, starts, index, element, numberOther);
4831	} else {
4832	//CoinPackedMatrix mm(matrix_);*
4833	matrix_->appendMinorFast(number, starts, index, element);
4834	//mm.appendRows(number,starts,index,element,numberOther);
4835	//if (!mm.isEquivalent(matrix_)) {*
4836	//printf("bad append\n");
4837	//abort();
4838	//}
4839	}
4840	} else {
4841	// columns
4842	if (!matrix_->isColOrdered() && numberOther > matrix_->getNumRows())
4843	matrix_->setDimensions(numberOther, -`1`);
4844	numberErrors = matrix_->appendCols(number, starts, index, element, numberOther);
4845	}
4846	clearCopies();
4847	numberActiveColumns_ = matrix_->getNumCols();
4848	return numberErrors;
4849	}
4850	void
4851	ClpPackedMatrix::specialRowCopy(ClpSimplex * model, const ClpMatrixBase * rowCopy)
4852	{
4853	delete rowCopy_;
4854	rowCopy_ = new ClpPackedMatrix2 (model, rowCopy->getPackedMatrix());
4855	// See if anything in it
4856	if (!rowCopy_->usefulInfo()) {
4857	delete rowCopy_;
4858	rowCopy_ = NULL;
4859	flags_ &= ~`4`;
4860	} else {
4861	flags_ \|= `4`;
4862	}
4863	}
4864	void
4865	ClpPackedMatrix::specialColumnCopy(ClpSimplex * model)
4866	{
4867	delete columnCopy_;
4868	if ((flags_ & `16`) != `0`) {
4869	columnCopy_ = new ClpPackedMatrix3 (model, matrix_);
4870	flags_ \|= `8`;
4871	} else {
4872	columnCopy_ = NULL;
4873	}
4874	}
4875	// Say we don't want special column copy
4876	void
4877	ClpPackedMatrix::releaseSpecialColumnCopy()
4878	{
4879	flags_ &= ~(`8` + `16`);
4880	delete columnCopy_;
4881	columnCopy_ = NULL;
4882	}
4883	// Correct sequence in and out to give true value
4884	void
4885	ClpPackedMatrix::correctSequence(const ClpSimplex * model, int & sequenceIn, int & sequenceOut)
4886	{
4887	if (columnCopy_) {
4888	if (sequenceIn != -`999`) {
4889	if (sequenceIn != sequenceOut) {
4890	if (sequenceIn < numberActiveColumns_)
4891	columnCopy_->swapOne(model, this, sequenceIn);
4892	if (sequenceOut < numberActiveColumns_)
4893	columnCopy_->swapOne(model, this, sequenceOut);
4894	}
4895	} else {
4896	// do all
4897	columnCopy_->sortBlocks(model);
4898	}
4899	}
4900	}
4901	// Check validity
4902	void
4903	ClpPackedMatrix::checkFlags(int type) const
4904	{
4905	int iColumn;
4906	// get matrix data pointers
4907	//const int row = matrix_->getIndices();*
4908	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
4909	const int * columnLength = matrix_->getVectorLengths();
4910	const double * elementByColumn = matrix_->getElements();
4911	if (!zeros()) {
4912	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
4913	CoinBigIndex j;
4914	for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) {
4915	if (!elementByColumn[j])
4916	abort();
4917	}
4918	}
4919	}
4920	if ((flags_ & `2`) == `0`) {
4921	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
4922	if (columnStart[iColumn+`1`] != columnStart[iColumn] + columnLength[iColumn]) {
4923	abort();
4924	}
4925	}
4926	}
4927	if (type) {
4928	if ((flags_ & `2`) != `0`) {
4929	bool ok = true;
4930	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
4931	if (columnStart[iColumn+`1`] != columnStart[iColumn] + columnLength[iColumn]) {
4932	ok = false;
4933	break;
4934	}
4935	}
4936	if (ok)
4937	COIN_DETAIL_PRINT(printf("flags_ could be 0\n"));
4938	}
4939	}
4940	}
4941	//#############################################################################
4942	// Constructors / Destructor / Assignment
4943	//#############################################################################
4944
4945	//-------------------------------------------------------------------
4946	// Default Constructor
4947	//-------------------------------------------------------------------
4948	ClpPackedMatrix2::ClpPackedMatrix2 ()
4949	: numberBlocks_(`0`),
4950	numberRows_(`0`),
4951	offset_(NULL),
4952	count_(NULL),
4953	rowStart_(NULL),
4954	column_(NULL),
4955	work_(NULL)
4956	{
4957	#ifdef THREAD
4958	threadId_ = NULL;
4959	info_ = NULL;
4960	#endif
4961	}
4962	//-------------------------------------------------------------------
4963	// Useful Constructor
4964	//-------------------------------------------------------------------
4965	ClpPackedMatrix2::ClpPackedMatrix2 (ClpSimplex * , const CoinPackedMatrix * rowCopy)
4966	: numberBlocks_(`0`),
4967	numberRows_(`0`),
4968	offset_(NULL),
4969	count_(NULL),
4970	rowStart_(NULL),
4971	column_(NULL),
4972	work_(NULL)
4973	{
4974	#ifdef THREAD
4975	threadId_ = NULL;
4976	info_ = NULL;
4977	#endif
4978	numberRows_ = rowCopy->getNumRows();
4979	if (!numberRows_)
4980	return;
4981	int numberColumns = rowCopy->getNumCols();
4982	const int * column = rowCopy->getIndices();
4983	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
4984	const int * length = rowCopy->getVectorLengths();
4985	const double * element = rowCopy->getElements();
4986	int chunk = `32768`; // tune
4987	//chunk=100;
4988	// tune
4989	#if 0
4990	int chunkY[`7`] = {`1024`, `2048`, `4096`, `8192`, `16384`, `32768`, `65535`};
4991	int its = model->maximumIterations();
4992	if (its >= `1000000` && its < `1000999`) {
4993	its -= `1000000`;
4994	its = its / `10`;
4995	if (its >= `7`) abort();
4996	chunk = chunkY[its];
4997	printf("chunk size %d\n", chunk);
4998	#define cpuid(func,ax,bx,cx,dx)\
4999	__asm__ __volatile__ ("cpuid":\
5000	"=a" (ax), "=b" (bx), "=c" (cx), "=d" (dx) : "a" (func));
5001	unsigned int a, b, c, d;
5002	int func = `0`;
5003	cpuid(func, a, b, c, d);
5004	{
5005	int i;
5006	unsigned int value;
5007	value = b;
5008	for (i = `0`; i < `4`; i++) {
5009	printf("%c", (value & `0xff`));
5010	value = value >> `8`;
5011	}
5012	value = d;
5013	for (i = `0`; i < `4`; i++) {
5014	printf("%c", (value & `0xff`));
5015	value = value >> `8`;
5016	}
5017	value = c;
5018	for (i = `0`; i < `4`; i++) {
5019	printf("%c", (value & `0xff`));
5020	value = value >> `8`;
5021	}
5022	printf("\n");
5023	int maxfunc = a;
5024	if (maxfunc > `10`) {
5025	printf("not intel?\n");
5026	abort();
5027	}
5028	for (func = `1`; func <= maxfunc; func++) {
5029	cpuid(func, a, b, c, d);
5030	printf("func %d, %x %x %x %x\n", func, a, b, c, d);
5031	}
5032	}
5033	#else
5034	if (numberColumns > `10000` \|\| chunk == `100`) {
5035	#endif
5036	} else {
5037	//printf("no chunk\n");
5038	return;
5039	}
5040	// Could also analyze matrix to get natural breaks
5041	numberBlocks_ = (numberColumns + chunk - `1`) / chunk;
5042	#ifdef THREAD
5043	// Get work areas
5044	threadId_ = new pthread_t [numberBlocks_];
5045	info_ = new dualColumn0Struct[numberBlocks_];
5046	#endif
5047	// Even out
5048	chunk = (numberColumns + numberBlocks_ - `1`) / numberBlocks_;
5049	offset_ = new int[numberBlocks_+`1`];
5050	offset_[numberBlocks_] = numberColumns;
5051	int nRow = numberBlocks_ * numberRows_;
5052	count_ = new unsigned short[nRow];
5053	memset(count_, `0`, nRow * sizeof(unsigned short));
5054	rowStart_ = new CoinBigIndex[nRow+numberRows_+`1`];
5055	CoinBigIndex nElement = rowStart[numberRows_];
5056	rowStart_[nRow+numberRows_] = nElement;
5057	column_ = new unsigned short [nElement];
5058	// assumes int <= double
5059	int sizeWork = `6` * numberBlocks_;
5060	work_ = new double[sizeWork];
5061	int iBlock;
5062	int nZero = `0`;
5063	for (iBlock = `0`; iBlock < numberBlocks_; iBlock++) {
5064	int start = iBlock * chunk;
5065	offset_[iBlock] = start;
5066	int end = start + chunk;
5067	for (int iRow = `0`; iRow < numberRows_; iRow++) {
5068	if (rowStart[iRow+`1`] != rowStart[iRow] + length[iRow]) {
5069	printf("not packed correctly - gaps\n");
5070	abort();
5071	}
5072	bool lastFound = false;
5073	int nFound = `0`;
5074	for (CoinBigIndex j = rowStart[iRow];
5075	j < rowStart[iRow] + length[iRow]; j++) {
5076	int iColumn = column[j];
5077	if (iColumn >= start) {
5078	if (iColumn < end) {
5079	if (!element[j]) {
5080	printf("not packed correctly - zero element\n");
5081	abort();
5082	}
5083	column_[j] = static_cast<unsigned short>(iColumn - start);
5084	nFound++;
5085	if (lastFound) {
5086	printf("not packed correctly - out of order\n");
5087	abort();
5088	}
5089	} else {
5090	//can't find any more
5091	lastFound = true;
5092	}
5093	}
5094	}
5095	count_[iRownumberBlocks_+iBlock] = static_cast<unsigned* short>(nFound);
5096	if (!nFound)
5097	nZero++;
5098	}
5099	}
5100	//double fraction = ((double) nZero)/((double) (numberBlocks_numberRows_));*
5101	//printf("%d empty blocks, %g%%\n",nZero,100.0fraction);*
5102	}
5103
5104	//-------------------------------------------------------------------
5105	// Copy constructor
5106	//-------------------------------------------------------------------
5107	ClpPackedMatrix2::ClpPackedMatrix2 (const ClpPackedMatrix2 & rhs)
5108	: numberBlocks_(rhs.numberBlocks_),
5109	numberRows_(rhs.numberRows_)
5110	{
5111	if (numberBlocks_) {
5112	offset_ = CoinCopyOfArray(rhs.offset_, numberBlocks_ + `1`);
5113	int nRow = numberBlocks_ * numberRows_;
5114	count_ = CoinCopyOfArray(rhs.count_, nRow);
5115	rowStart_ = CoinCopyOfArray(rhs.rowStart_, nRow + numberRows_ + `1`);
5116	CoinBigIndex nElement = rowStart_[nRow+numberRows_];
5117	column_ = CoinCopyOfArray(rhs.column_, nElement);
5118	int sizeWork = `6` * numberBlocks_;
5119	work_ = CoinCopyOfArray(rhs.work_, sizeWork);
5120	#ifdef THREAD
5121	threadId_ = new pthread_t [numberBlocks_];
5122	info_ = new dualColumn0Struct[numberBlocks_];
5123	#endif
5124	} else {
5125	offset_ = NULL;
5126	count_ = NULL;
5127	rowStart_ = NULL;
5128	column_ = NULL;
5129	work_ = NULL;
5130	#ifdef THREAD
5131	threadId_ = NULL;
5132	info_ = NULL;
5133	#endif
5134	}
5135	}
5136	//-------------------------------------------------------------------
5137	// Destructor
5138	//-------------------------------------------------------------------
5139	ClpPackedMatrix2::~ClpPackedMatrix2 ()
5140	{
5141	delete [] offset_;
5142	delete [] count_;
5143	delete [] rowStart_;
5144	delete [] column_;
5145	delete [] work_;
5146	#ifdef THREAD
5147	delete [] threadId_;
5148	delete [] info_;
5149	#endif
5150	}
5151
5152	//----------------------------------------------------------------
5153	// Assignment operator
5154	//-------------------------------------------------------------------
5155	ClpPackedMatrix2 &
5156	ClpPackedMatrix2::operator=(const ClpPackedMatrix2& rhs)
5157	{
5158	if (this != &rhs) {
5159	numberBlocks_ = rhs.numberBlocks_;
5160	numberRows_ = rhs.numberRows_;
5161	delete [] offset_;
5162	delete [] count_;
5163	delete [] rowStart_;
5164	delete [] column_;
5165	delete [] work_;
5166	#ifdef THREAD
5167	delete [] threadId_;
5168	delete [] info_;
5169	#endif
5170	if (numberBlocks_) {
5171	offset_ = CoinCopyOfArray(rhs.offset_, numberBlocks_ + `1`);
5172	int nRow = numberBlocks_ * numberRows_;
5173	count_ = CoinCopyOfArray(rhs.count_, nRow);
5174	rowStart_ = CoinCopyOfArray(rhs.rowStart_, nRow + numberRows_ + `1`);
5175	CoinBigIndex nElement = rowStart_[nRow+numberRows_];
5176	column_ = CoinCopyOfArray(rhs.column_, nElement);
5177	int sizeWork = `6` * numberBlocks_;
5178	work_ = CoinCopyOfArray(rhs.work_, sizeWork);
5179	#ifdef THREAD
5180	threadId_ = new pthread_t [numberBlocks_];
5181	info_ = new dualColumn0Struct[numberBlocks_];
5182	#endif
5183	} else {
5184	offset_ = NULL;
5185	count_ = NULL;
5186	rowStart_ = NULL;
5187	column_ = NULL;
5188	work_ = NULL;
5189	#ifdef THREAD
5190	threadId_ = NULL;
5191	info_ = NULL;
5192	#endif
5193	}
5194	}
5195	return *this;
5196	}
5197	static int dualColumn0(const ClpSimplex * model, double * spare,
5198	int * spareIndex, const double * arrayTemp,
5199	const int * indexTemp, int numberIn,
5200	int offset, double acceptablePivot, double * bestPossiblePtr,
5201	double * upperThetaPtr, int * posFreePtr, double * freePivotPtr)
5202	{
5203	// do dualColumn0
5204	int i;
5205	int numberRemaining = `0`;
5206	double bestPossible = `0.0`;
5207	double upperTheta = `1.0e31`;
5208	double freePivot = acceptablePivot;
5209	int posFree = -`1`;
5210	const double * reducedCost = model->djRegion(`1`);
5211	double dualTolerance = model->dualTolerance();
5212	// We can also see if infeasible or pivoting on free
5213	double tentativeTheta = `1.0e25`;
5214	for (i = `0`; i < numberIn; i++) {
5215	double alpha = arrayTemp[i];
5216	int iSequence = indexTemp[i] + offset;
5217	double oldValue;
5218	double value;
5219	bool keep;
5220
5221	switch(model->getStatus(iSequence)) {
5222
5223	case ClpSimplex::basic:
5224	case ClpSimplex::isFixed:
5225	break;
5226	case ClpSimplex::isFree:
5227	case ClpSimplex::superBasic:
5228	bestPossible = CoinMax(bestPossible, fabs(alpha));
5229	oldValue = reducedCost[iSequence];
5230	// If free has to be very large - should come in via dualRow
5231	if (model->getStatus(iSequence) == ClpSimplex::isFree && fabs(alpha) < `1.0e-3`)
5232	break;
5233	if (oldValue > dualTolerance) {
5234	keep = true;
5235	} else if (oldValue < -dualTolerance) {
5236	keep = true;
5237	} else {
5238	if (fabs(alpha) > CoinMax(`10.0` * acceptablePivot, `1.0e-5`))
5239	keep = true;
5240	else
5241	keep = false;
5242	}
5243	if (keep) {
5244	// free - choose largest
5245	if (fabs(alpha) > freePivot) {
5246	freePivot = fabs(alpha);
5247	posFree = i;
5248	}
5249	}
5250	break;
5251	case ClpSimplex::atUpperBound:
5252	oldValue = reducedCost[iSequence];
5253	value = oldValue - tentativeTheta * alpha;
5254	//assert (oldValue<=dualTolerance1.0001);*
5255	if (value > dualTolerance) {
5256	bestPossible = CoinMax(bestPossible, -alpha);
5257	value = oldValue - upperTheta * alpha;
5258	if (value > dualTolerance && -alpha >= acceptablePivot)
5259	upperTheta = (oldValue - dualTolerance) / alpha;
5260	// add to list
5261	spare[numberRemaining] = alpha;
5262	spareIndex[numberRemaining++] = iSequence;
5263	}
5264	break;
5265	case ClpSimplex::atLowerBound:
5266	oldValue = reducedCost[iSequence];
5267	value = oldValue - tentativeTheta * alpha;
5268	//assert (oldValue>=-dualTolerance1.0001);*
5269	if (value < -dualTolerance) {
5270	bestPossible = CoinMax(bestPossible, alpha);
5271	value = oldValue - upperTheta * alpha;
5272	if (value < -dualTolerance && alpha >= acceptablePivot)
5273	upperTheta = (oldValue + dualTolerance) / alpha;
5274	// add to list
5275	spare[numberRemaining] = alpha;
5276	spareIndex[numberRemaining++] = iSequence;
5277	}
5278	break;
5279	}
5280	}
5281	*bestPossiblePtr = bestPossible;
5282	*upperThetaPtr = upperTheta;
5283	*freePivotPtr = freePivot;
5284	*posFreePtr = posFree;
5285	return numberRemaining;
5286	}
5287	static int doOneBlock(double * array, int * index,
5288	const double * pi, const CoinBigIndex * rowStart, const double * element,
5289	const unsigned short * column, int numberInRowArray, int numberLook)
5290	{
5291	int iWhich = `0`;
5292	int nextN = `0`;
5293	CoinBigIndex nextStart = `0`;
5294	double nextPi = `0.0`;
5295	for (; iWhich < numberInRowArray; iWhich++) {
5296	nextStart = rowStart[`0`];
5297	nextN = rowStart[numberInRowArray] - nextStart;
5298	rowStart++;
5299	if (nextN) {
5300	nextPi = pi[iWhich];
5301	break;
5302	}
5303	}
5304	int i;
5305	while (iWhich < numberInRowArray) {
5306	double value = nextPi;
5307	CoinBigIndex j = nextStart;
5308	int n = nextN;
5309	// get next
5310	iWhich++;
5311	for (; iWhich < numberInRowArray; iWhich++) {
5312	nextStart = rowStart[`0`];
5313	nextN = rowStart[numberInRowArray] - nextStart;
5314	rowStart++;
5315	if (nextN) {
5316	//coin_prefetch_const(element + nextStart);
5317	nextPi = pi[iWhich];
5318	break;
5319	}
5320	}
5321	CoinBigIndex end = j + n;
5322	//coin_prefetch_const(element+rowStart_[i+1]);
5323	//coin_prefetch_const(column_+rowStart_[i+1]);
5324	if (n < `100`) {
5325	if ((n & `1`) != `0`) {
5326	unsigned int jColumn = column[j];
5327	array[jColumn] -= value * element[j];
5328	j++;
5329	}
5330	//coin_prefetch_const(column + nextStart);
5331	for (; j < end; j += `2`) {
5332	unsigned int jColumn0 = column[j];
5333	double value0 = value * element[j];
5334	unsigned int jColumn1 = column[j+`1`];
5335	double value1 = value * element[j+`1`];
5336	array[jColumn0] -= value0;
5337	array[jColumn1] -= value1;
5338	}
5339	} else {
5340	if ((n & `1`) != `0`) {
5341	unsigned int jColumn = column[j];
5342	array[jColumn] -= value * element[j];
5343	j++;
5344	}
5345	if ((n & `2`) != `0`) {
5346	unsigned int jColumn0 = column[j];
5347	double value0 = value * element[j];
5348	unsigned int jColumn1 = column[j+`1`];
5349	double value1 = value * element[j+`1`];
5350	array[jColumn0] -= value0;
5351	array[jColumn1] -= value1;
5352	j += `2`;
5353	}
5354	if ((n & `4`) != `0`) {
5355	unsigned int jColumn0 = column[j];
5356	double value0 = value * element[j];
5357	unsigned int jColumn1 = column[j+`1`];
5358	double value1 = value * element[j+`1`];
5359	unsigned int jColumn2 = column[j+`2`];
5360	double value2 = value * element[j+`2`];
5361	unsigned int jColumn3 = column[j+`3`];
5362	double value3 = value * element[j+`3`];
5363	array[jColumn0] -= value0;
5364	array[jColumn1] -= value1;
5365	array[jColumn2] -= value2;
5366	array[jColumn3] -= value3;
5367	j += `4`;
5368	}
5369	//coin_prefetch_const(column+nextStart);
5370	for (; j < end; j += `8`) {
5371	//coin_prefetch_const(element + j + 16);
5372	unsigned int jColumn0 = column[j];
5373	double value0 = value * element[j];
5374	unsigned int jColumn1 = column[j+`1`];
5375	double value1 = value * element[j+`1`];
5376	unsigned int jColumn2 = column[j+`2`];
5377	double value2 = value * element[j+`2`];
5378	unsigned int jColumn3 = column[j+`3`];
5379	double value3 = value * element[j+`3`];
5380	array[jColumn0] -= value0;
5381	array[jColumn1] -= value1;
5382	array[jColumn2] -= value2;
5383	array[jColumn3] -= value3;
5384	//coin_prefetch_const(column + j + 16);
5385	jColumn0 = column[j+`4`];
5386	value0 = value * element[j+`4`];
5387	jColumn1 = column[j+`5`];
5388	value1 = value * element[j+`5`];
5389	jColumn2 = column[j+`6`];
5390	value2 = value * element[j+`6`];
5391	jColumn3 = column[j+`7`];
5392	value3 = value * element[j+`7`];
5393	array[jColumn0] -= value0;
5394	array[jColumn1] -= value1;
5395	array[jColumn2] -= value2;
5396	array[jColumn3] -= value3;
5397	}
5398	}
5399	}
5400	// get rid of tiny values
5401	int nSmall = numberLook;
5402	int numberNonZero = `0`;
5403	for (i = `0`; i < nSmall; i++) {
5404	double value = array[i];
5405	array[i] = `0.0`;
5406	if (fabs(value) > `1.0e-12`) {
5407	array[numberNonZero] = value;
5408	index[numberNonZero++] = i;
5409	}
5410	}
5411	for (; i < numberLook; i += `4`) {
5412	double value0 = array[i+`0`];
5413	double value1 = array[i+`1`];
5414	double value2 = array[i+`2`];
5415	double value3 = array[i+`3`];
5416	array[i+`0`] = `0.0`;
5417	array[i+`1`] = `0.0`;
5418	array[i+`2`] = `0.0`;
5419	array[i+`3`] = `0.0`;
5420	if (fabs(value0) > `1.0e-12`) {
5421	array[numberNonZero] = value0;
5422	index[numberNonZero++] = i + `0`;
5423	}
5424	if (fabs(value1) > `1.0e-12`) {
5425	array[numberNonZero] = value1;
5426	index[numberNonZero++] = i + `1`;
5427	}
5428	if (fabs(value2) > `1.0e-12`) {
5429	array[numberNonZero] = value2;
5430	index[numberNonZero++] = i + `2`;
5431	}
5432	if (fabs(value3) > `1.0e-12`) {
5433	array[numberNonZero] = value3;
5434	index[numberNonZero++] = i + `3`;
5435	}
5436	}
5437	return numberNonZero;
5438	}
5439	#ifdef THREAD
5440	static void * doOneBlockThread(void * voidInfo)
5441	{
5442	dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
5443	*(info->numberInPtr) = doOneBlock(info->arrayTemp, info->indexTemp, info->pi,
5444	info->rowStart, info->element, info->column,
5445	info->numberInRowArray, info->numberLook);
5446	return NULL;
5447	}
5448	static void * doOneBlockAnd0Thread(void * voidInfo)
5449	{
5450	dualColumn0Struct * info = (dualColumn0Struct *) voidInfo;
5451	*(info->numberInPtr) = doOneBlock(info->arrayTemp, info->indexTemp, info->pi,
5452	info->rowStart, info->element, info->column,
5453	info->numberInRowArray, info->numberLook);
5454	*(info->numberOutPtr) = dualColumn0(info->model, info->spare,
5455	info->spareIndex, (const double *)info->arrayTemp,
5456	(const int ) info->indexTemp, (info->numberInPtr),
5457	info->offset, info->acceptablePivot, info->bestPossiblePtr,
5458	info->upperThetaPtr, info->posFreePtr, info->freePivotPtr);
5459	return NULL;
5460	}
5461	#endif
5462	/ Return <code>x * scalar * A in <code>z</code>.*
5463	Note - x packed and z will be packed mode
5464	Squashes small elements and knows about ClpSimplex /*
5465	void
5466	ClpPackedMatrix2::transposeTimes(const ClpSimplex * model,
5467	const CoinPackedMatrix * rowCopy,
5468	const CoinIndexedVector * rowArray,
5469	CoinIndexedVector * spareArray,
5470	CoinIndexedVector * columnArray) const
5471	{
5472	// See if dualColumn0 coding wanted
5473	bool dualColumn = model->spareIntArray_[`0`] == `1`;
5474	double acceptablePivot = model->spareDoubleArray_[`0`];
5475	double bestPossible = `0.0`;
5476	double upperTheta = `1.0e31`;
5477	double freePivot = acceptablePivot;
5478	int posFree = -`1`;
5479	int numberRemaining = `0`;
5480	//if (model->numberIterations()>=200000) {
5481	//printf("time %g\n",CoinCpuTime()-startTime);
5482	//exit(77);
5483	//}
5484	double * pi = rowArray->denseVector();
5485	int numberNonZero = `0`;
5486	int * index = columnArray->getIndices();
5487	double * array = columnArray->denseVector();
5488	int numberInRowArray = rowArray->getNumElements();
5489	const int * whichRow = rowArray->getIndices();
5490	double * element = const_cast<double *>(rowCopy->getElements());
5491	const CoinBigIndex * rowStart = rowCopy->getVectorStarts();
5492	int i;
5493	CoinBigIndex * rowStart2 = rowStart_;
5494	if (!dualColumn) {
5495	for (i = `0`; i < numberInRowArray; i++) {
5496	int iRow = whichRow[i];
5497	CoinBigIndex start = rowStart[iRow];
5498	*rowStart2 = start;
5499	unsigned short * count1 = count_ + iRow * numberBlocks_;
5500	int put = `0`;
5501	for (int j = `0`; j < numberBlocks_; j++) {
5502	put += numberInRowArray;
5503	start += count1[j];
5504	rowStart2[put] = start;
5505	}
5506	rowStart2 ++;
5507	}
5508	} else {
5509	// also do dualColumn stuff
5510	double * spare = spareArray->denseVector();
5511	int * spareIndex = spareArray->getIndices();
5512	const double * reducedCost = model->djRegion(`0`);
5513	double dualTolerance = model->dualTolerance();
5514	// We can also see if infeasible or pivoting on free
5515	double tentativeTheta = `1.0e25`;
5516	int addSequence = model->numberColumns();
5517	for (i = `0`; i < numberInRowArray; i++) {
5518	int iRow = whichRow[i];
5519	double alpha = pi[i];
5520	double oldValue;
5521	double value;
5522	bool keep;
5523
5524	switch(model->getStatus(iRow + addSequence)) {
5525
5526	case ClpSimplex::basic:
5527	case ClpSimplex::isFixed:
5528	break;
5529	case ClpSimplex::isFree:
5530	case ClpSimplex::superBasic:
5531	bestPossible = CoinMax(bestPossible, fabs(alpha));
5532	oldValue = reducedCost[iRow];
5533	// If free has to be very large - should come in via dualRow
5534	if (model->getStatus(iRow + addSequence) == ClpSimplex::isFree && fabs(alpha) < `1.0e-3`)
5535	break;
5536	if (oldValue > dualTolerance) {
5537	keep = true;
5538	} else if (oldValue < -dualTolerance) {
5539	keep = true;
5540	} else {
5541	if (fabs(alpha) > CoinMax(`10.0` * acceptablePivot, `1.0e-5`))
5542	keep = true;
5543	else
5544	keep = false;
5545	}
5546	if (keep) {
5547	// free - choose largest
5548	if (fabs(alpha) > freePivot) {
5549	freePivot = fabs(alpha);
5550	posFree = i + addSequence;
5551	}
5552	}
5553	break;
5554	case ClpSimplex::atUpperBound:
5555	oldValue = reducedCost[iRow];
5556	value = oldValue - tentativeTheta * alpha;
5557	//assert (oldValue<=dualTolerance1.0001);*
5558	if (value > dualTolerance) {
5559	bestPossible = CoinMax(bestPossible, -alpha);
5560	value = oldValue - upperTheta * alpha;
5561	if (value > dualTolerance && -alpha >= acceptablePivot)
5562	upperTheta = (oldValue - dualTolerance) / alpha;
5563	// add to list
5564	spare[numberRemaining] = alpha;
5565	spareIndex[numberRemaining++] = iRow + addSequence;
5566	}
5567	break;
5568	case ClpSimplex::atLowerBound:
5569	oldValue = reducedCost[iRow];
5570	value = oldValue - tentativeTheta * alpha;
5571	//assert (oldValue>=-dualTolerance1.0001);*
5572	if (value < -dualTolerance) {
5573	bestPossible = CoinMax(bestPossible, alpha);
5574	value = oldValue - upperTheta * alpha;
5575	if (value < -dualTolerance && alpha >= acceptablePivot)
5576	upperTheta = (oldValue + dualTolerance) / alpha;
5577	// add to list
5578	spare[numberRemaining] = alpha;
5579	spareIndex[numberRemaining++] = iRow + addSequence;
5580	}
5581	break;
5582	}
5583	CoinBigIndex start = rowStart[iRow];
5584	*rowStart2 = start;
5585	unsigned short * count1 = count_ + iRow * numberBlocks_;
5586	int put = `0`;
5587	for (int j = `0`; j < numberBlocks_; j++) {
5588	put += numberInRowArray;
5589	start += count1[j];
5590	rowStart2[put] = start;
5591	}
5592	rowStart2 ++;
5593	}
5594	}
5595	double * spare = spareArray->denseVector();
5596	int * spareIndex = spareArray->getIndices();
5597	int saveNumberRemaining = numberRemaining;
5598	int iBlock;
5599	for (iBlock = `0`; iBlock < numberBlocks_; iBlock++) {
5600	double * dwork = work_ + `6` * iBlock;
5601	int * iwork = reinterpret_cast<int *> (dwork + `3`);
5602	if (!dualColumn) {
5603	#ifndef THREAD
5604	int offset = offset_[iBlock];
5605	int offset3 = offset;
5606	offset = numberNonZero;
5607	double * arrayTemp = array + offset;
5608	int * indexTemp = index + offset;
5609	iwork[`0`] = doOneBlock(arrayTemp, indexTemp, pi, rowStart_ + numberInRowArray * iBlock,
5610	element, column_, numberInRowArray, offset_[iBlock+`1`] - offset);
5611	int number = iwork[`0`];
5612	for (i = `0`; i < number; i++) {
5613	//double value = arrayTemp[i];
5614	//arrayTemp[i]=0.0;
5615	//array[numberNonZero]=value;
5616	index[numberNonZero++] = indexTemp[i] + offset3;
5617	}
5618	#else
5619	int offset = offset_[iBlock];
5620	double * arrayTemp = array + offset;
5621	int * indexTemp = index + offset;
5622	dualColumn0Struct * infoPtr = info_ + iBlock;
5623	infoPtr->arrayTemp = arrayTemp;
5624	infoPtr->indexTemp = indexTemp;
5625	infoPtr->numberInPtr = &iwork[`0`];
5626	infoPtr->pi = pi;
5627	infoPtr->rowStart = rowStart_ + numberInRowArray * iBlock;
5628	infoPtr->element = element;
5629	infoPtr->column = column_;
5630	infoPtr->numberInRowArray = numberInRowArray;
5631	infoPtr->numberLook = offset_[iBlock+`1`] - offset;
5632	pthread_create(&threadId_[iBlock], NULL, doOneBlockThread, infoPtr);
5633	#endif
5634	} else {
5635	#ifndef THREAD
5636	int offset = offset_[iBlock];
5637	// allow for already saved
5638	int offset2 = offset + saveNumberRemaining;
5639	int offset3 = offset;
5640	offset = numberNonZero;
5641	offset2 = numberRemaining;
5642	double * arrayTemp = array + offset;
5643	int * indexTemp = index + offset;
5644	iwork[`0`] = doOneBlock(arrayTemp, indexTemp, pi, rowStart_ + numberInRowArray * iBlock,
5645	element, column_, numberInRowArray, offset_[iBlock+`1`] - offset);
5646	iwork[`1`] = dualColumn0(model, spare + offset2,
5647	spareIndex + offset2,
5648	arrayTemp, indexTemp,
5649	iwork[`0`], offset3, acceptablePivot,
5650	&dwork[`0`], &dwork[`1`], &iwork[`2`],
5651	&dwork[`2`]);
5652	int number = iwork[`0`];
5653	int numberLook = iwork[`1`];
5654	#if 1
5655	numberRemaining += numberLook;
5656	#else
5657	double * spareTemp = spare + offset2;
5658	const int * spareIndexTemp = spareIndex + offset2;
5659	for (i = `0`; i < numberLook; i++) {
5660	double value = spareTemp[i];
5661	spareTemp[i] = `0.0`;
5662	spare[numberRemaining] = value;
5663	spareIndex[numberRemaining++] = spareIndexTemp[i];
5664	}
5665	#endif
5666	if (dwork[`2`] > freePivot) {
5667	freePivot = dwork[`2`];
5668	posFree = iwork[`2`] + numberNonZero;
5669	}
5670	upperTheta = CoinMin(dwork[`1`], upperTheta);
5671	bestPossible = CoinMax(dwork[`0`], bestPossible);
5672	for (i = `0`; i < number; i++) {
5673	// double value = arrayTemp[i];
5674	//arrayTemp[i]=0.0;
5675	//array[numberNonZero]=value;
5676	index[numberNonZero++] = indexTemp[i] + offset3;
5677	}
5678	#else
5679	int offset = offset_[iBlock];
5680	// allow for already saved
5681	int offset2 = offset + saveNumberRemaining;
5682	double * arrayTemp = array + offset;
5683	int * indexTemp = index + offset;
5684	dualColumn0Struct * infoPtr = info_ + iBlock;
5685	infoPtr->model = model;
5686	infoPtr->spare = spare + offset2;
5687	infoPtr->spareIndex = spareIndex + offset2;
5688	infoPtr->arrayTemp = arrayTemp;
5689	infoPtr->indexTemp = indexTemp;
5690	infoPtr->numberInPtr = &iwork[`0`];
5691	infoPtr->offset = offset;
5692	infoPtr->acceptablePivot = acceptablePivot;
5693	infoPtr->bestPossiblePtr = &dwork[`0`];
5694	infoPtr->upperThetaPtr = &dwork[`1`];
5695	infoPtr->posFreePtr = &iwork[`2`];
5696	infoPtr->freePivotPtr = &dwork[`2`];
5697	infoPtr->numberOutPtr = &iwork[`1`];
5698	infoPtr->pi = pi;
5699	infoPtr->rowStart = rowStart_ + numberInRowArray * iBlock;
5700	infoPtr->element = element;
5701	infoPtr->column = column_;
5702	infoPtr->numberInRowArray = numberInRowArray;
5703	infoPtr->numberLook = offset_[iBlock+`1`] - offset;
5704	if (iBlock >= `2`)
5705	pthread_join(threadId_[iBlock-`2`], NULL);
5706	pthread_create(threadId_ + iBlock, NULL, doOneBlockAnd0Thread, infoPtr);
5707	//pthread_join(threadId_[iBlock],NULL);
5708	#endif
5709	}
5710	}
5711	for ( iBlock = CoinMax(`0`, numberBlocks_ - `2`); iBlock < numberBlocks_; iBlock++) {
5712	#ifdef THREAD
5713	pthread_join(threadId_[iBlock], NULL);
5714	#endif
5715	}
5716	#ifdef THREAD
5717	for ( iBlock = `0`; iBlock < numberBlocks_; iBlock++) {
5718	//pthread_join(threadId_[iBlock],NULL);
5719	int offset = offset_[iBlock];
5720	double * dwork = work_ + `6` * iBlock;
5721	int * iwork = (int *) (dwork + `3`);
5722	int number = iwork[`0`];
5723	if (dualColumn) {
5724	// allow for already saved
5725	int offset2 = offset + saveNumberRemaining;
5726	int numberLook = iwork[`1`];
5727	double * spareTemp = spare + offset2;
5728	const int * spareIndexTemp = spareIndex + offset2;
5729	for (i = `0`; i < numberLook; i++) {
5730	double value = spareTemp[i];
5731	spareTemp[i] = `0.0`;
5732	spare[numberRemaining] = value;
5733	spareIndex[numberRemaining++] = spareIndexTemp[i];
5734	}
5735	if (dwork[`2`] > freePivot) {
5736	freePivot = dwork[`2`];
5737	posFree = iwork[`2`] + numberNonZero;
5738	}
5739	upperTheta = CoinMin(dwork[`1`], upperTheta);
5740	bestPossible = CoinMax(dwork[`0`], bestPossible);
5741	}
5742	double * arrayTemp = array + offset;
5743	const int * indexTemp = index + offset;
5744	for (i = `0`; i < number; i++) {
5745	double value = arrayTemp[i];
5746	arrayTemp[i] = `0.0`;
5747	array[numberNonZero] = value;
5748	index[numberNonZero++] = indexTemp[i] + offset;
5749	}
5750	}
5751	#endif
5752	columnArray->setNumElements(numberNonZero);
5753	columnArray->setPackedMode(true);
5754	if (dualColumn) {
5755	model->spareDoubleArray_[`0`] = upperTheta;
5756	model->spareDoubleArray_[`1`] = bestPossible;
5757	// and theta and alpha and sequence
5758	if (posFree < `0`) {
5759	model->spareIntArray_[`1`] = -`1`;
5760	} else {
5761	const double * reducedCost = model->djRegion(`0`);
5762	double alpha;
5763	int numberColumns = model->numberColumns();
5764	if (posFree < numberColumns) {
5765	alpha = columnArray->denseVector()[posFree];
5766	posFree = columnArray->getIndices()[posFree];
5767	} else {
5768	alpha = rowArray->denseVector()[posFree-numberColumns];
5769	posFree = rowArray->getIndices()[posFree-numberColumns] + numberColumns;
5770	}
5771	model->spareDoubleArray_[`2`] = fabs(reducedCost[posFree] / alpha);
5772	model->spareDoubleArray_[`3`] = alpha;
5773	model->spareIntArray_[`1`] = posFree;
5774	}
5775	spareArray->setNumElements(numberRemaining);
5776	// signal done
5777	model->spareIntArray_[`0`] = -`1`;
5778	}
5779	}
5780	/ Default constructor. /
5781	ClpPackedMatrix3::ClpPackedMatrix3()
5782	: numberBlocks_(`0`),
5783	numberColumns_(`0`),
5784	column_(NULL),
5785	start_(NULL),
5786	row_(NULL),
5787	element_(NULL),
5788	block_(NULL)
5789	{
5790	}
5791	/ Constructor from copy. /
5792	ClpPackedMatrix3::ClpPackedMatrix3(ClpSimplex * model, const CoinPackedMatrix * columnCopy)
5793	: numberBlocks_(`0`),
5794	numberColumns_(`0`),
5795	column_(NULL),
5796	start_(NULL),
5797	row_(NULL),
5798	element_(NULL),
5799	block_(NULL)
5800	{
5801	#define MINBLOCK 6
5802	#define MAXBLOCK 100
5803	#define MAXUNROLL 10
5804	numberColumns_ = model->getNumCols();
5805	int numberColumns = columnCopy->getNumCols();
5806	assert (numberColumns_ >= numberColumns);
5807	int numberRows = columnCopy->getNumRows();
5808	int * counts = new int[numberRows+`1`];
5809	CoinZeroN(counts, numberRows + `1`);
5810	CoinBigIndex nels = `0`;
5811	CoinBigIndex nZeroEl = `0`;
5812	int iColumn;
5813	// get matrix data pointers
5814	const int * row = columnCopy->getIndices();
5815	const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
5816	const int * columnLength = columnCopy->getVectorLengths();
5817	const double * elementByColumn = columnCopy->getElements();
5818	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
5819	CoinBigIndex start = columnStart[iColumn];
5820	int n = columnLength[iColumn];
5821	CoinBigIndex end = start + n;
5822	int kZero = `0`;
5823	for (CoinBigIndex j = start; j < end; j++) {
5824	if(!elementByColumn[j])
5825	kZero++;
5826	}
5827	n -= kZero;
5828	nZeroEl += kZero;
5829	nels += n;
5830	counts[n]++;
5831	}
5832	counts[`0`] += numberColumns_ - numberColumns;
5833	int nZeroColumns = counts[`0`];
5834	counts[`0`] = -`1`;
5835	numberColumns_ -= nZeroColumns;
5836	column_ = new int [`2`*numberColumns_+nZeroColumns];
5837	int * lookup = column_ + numberColumns_;
5838	row_ = new int[nels];
5839	element_ = new double[nels];
5840	int nOdd = `0`;
5841	CoinBigIndex nInOdd = `0`;
5842	int i;
5843	for (i = `1`; i <= numberRows; i++) {
5844	int n = counts[i];
5845	if (n) {
5846	if (n < MINBLOCK \|\| i > MAXBLOCK) {
5847	nOdd += n;
5848	counts[i] = -`1`;
5849	nInOdd += n * i;
5850	} else {
5851	numberBlocks_++;
5852	}
5853	} else {
5854	counts[i] = -`1`;
5855	}
5856	}
5857	start_ = new CoinBigIndex [nOdd+`1`];
5858	// even if no blocks do a dummy one
5859	numberBlocks_ = CoinMax(numberBlocks_, `1`);
5860	block_ = new blockStruct [numberBlocks_];
5861	memset(block_, `0`, numberBlocks_ * sizeof(blockStruct));
5862	// Fill in what we can
5863	int nTotal = nOdd;
5864	// in case no blocks
5865	block_->startIndices_ = nTotal;
5866	nels = nInOdd;
5867	int nBlock = `0`;
5868	for (i = `0`; i <= CoinMin(MAXBLOCK, numberRows); i++) {
5869	if (counts[i] > `0`) {
5870	blockStruct * block = block_ + nBlock;
5871	int n = counts[i];
5872	counts[i] = nBlock; // backward pointer
5873	nBlock++;
5874	block->startIndices_ = nTotal;
5875	block->startElements_ = nels;
5876	block->numberElements_ = i;
5877	// up counts
5878	nTotal += n;
5879	nels += n * i;
5880	}
5881	}
5882	for (iColumn = numberColumns; iColumn < numberColumns_; iColumn++)
5883	lookup[iColumn] = -`1`;
5884	// fill
5885	start_[`0`] = `0`;
5886	nOdd = `0`;
5887	nInOdd = `0`;
5888	const double * columnScale = model->columnScale();
5889	for (iColumn = `0`; iColumn < numberColumns; iColumn++) {
5890	CoinBigIndex start = columnStart[iColumn];
5891	int n = columnLength[iColumn];
5892	CoinBigIndex end = start + n;
5893	int kZero = `0`;
5894	for (CoinBigIndex j = start; j < end; j++) {
5895	if(!elementByColumn[j])
5896	kZero++;
5897	}
5898	n -= kZero;
5899	if (n) {
5900	int iBlock = counts[n];
5901	if (iBlock >= `0`) {
5902	blockStruct * block = block_ + iBlock;
5903	int k = block->numberInBlock_;
5904	block->numberInBlock_ ++;
5905	column_[block->startIndices_+k] = iColumn;
5906	lookup[iColumn] = k;
5907	CoinBigIndex put = block->startElements_ + k * n;
5908	for (CoinBigIndex j = start; j < end; j++) {
5909	double value = elementByColumn[j];
5910	if(value) {
5911	if (columnScale)
5912	value *= columnScale[iColumn];
5913	element_[put] = value;
5914	row_[put++] = row[j];
5915	}
5916	}
5917	} else {
5918	// odd ones
5919	for (CoinBigIndex j = start; j < end; j++) {
5920	double value = elementByColumn[j];
5921	if(value) {
5922	if (columnScale)
5923	value *= columnScale[iColumn];
5924	element_[nInOdd] = value;
5925	row_[nInOdd++] = row[j];
5926	}
5927	}
5928	column_[nOdd] = iColumn;
5929	lookup[iColumn] = -`1`;
5930	nOdd++;
5931	start_[nOdd] = nInOdd;
5932	}
5933	} else {
5934	// zero column
5935	lookup[iColumn] = -`1`;
5936	}
5937	}
5938	delete [] counts;
5939	}
5940	/ Destructor /
5941	ClpPackedMatrix3::~ClpPackedMatrix3()
5942	{
5943	delete [] column_;
5944	delete [] start_;
5945	delete [] row_;
5946	delete [] element_;
5947	delete [] block_;
5948	}
5949	/ The copy constructor. /
5950	ClpPackedMatrix3::ClpPackedMatrix3(const ClpPackedMatrix3 & rhs)
5951	: numberBlocks_(rhs.numberBlocks_),
5952	numberColumns_(rhs.numberColumns_),
5953	column_(NULL),
5954	start_(NULL),
5955	row_(NULL),
5956	element_(NULL),
5957	block_(NULL)
5958	{
5959	if (rhs.numberBlocks_) {
5960	block_ = CoinCopyOfArray(rhs.block_, numberBlocks_);
5961	column_ = CoinCopyOfArray(rhs.column_, `2` * numberColumns_);
5962	int numberOdd = block_->startIndices_;
5963	start_ = CoinCopyOfArray(rhs.start_, numberOdd + `1`);
5964	blockStruct * lastBlock = block_ + (numberBlocks_ - `1`);
5965	CoinBigIndex numberElements = lastBlock->startElements_ +
5966	lastBlock->numberInBlock_ * lastBlock->numberElements_;
5967	row_ = CoinCopyOfArray(rhs.row_, numberElements);
5968	element_ = CoinCopyOfArray(rhs.element_, numberElements);
5969	}
5970	}
5971	ClpPackedMatrix3&
5972	ClpPackedMatrix3::operator=(const ClpPackedMatrix3 & rhs)
5973	{
5974	if (this != &rhs) {
5975	delete [] column_;
5976	delete [] start_;
5977	delete [] row_;
5978	delete [] element_;
5979	delete [] block_;
5980	numberBlocks_ = rhs.numberBlocks_;
5981	numberColumns_ = rhs.numberColumns_;
5982	if (rhs.numberBlocks_) {
5983	block_ = CoinCopyOfArray(rhs.block_, numberBlocks_);
5984	column_ = CoinCopyOfArray(rhs.column_, `2` * numberColumns_);
5985	int numberOdd = block_->startIndices_;
5986	start_ = CoinCopyOfArray(rhs.start_, numberOdd + `1`);
5987	blockStruct * lastBlock = block_ + (numberBlocks_ - `1`);
5988	CoinBigIndex numberElements = lastBlock->startElements_ +
5989	lastBlock->numberInBlock_ * lastBlock->numberElements_;
5990	row_ = CoinCopyOfArray(rhs.row_, numberElements);
5991	element_ = CoinCopyOfArray(rhs.element_, numberElements);
5992	} else {
5993	column_ = NULL;
5994	start_ = NULL;
5995	row_ = NULL;
5996	element_ = NULL;
5997	block_ = NULL;
5998	}
5999	}
6000	return *this;
6001	}
6002	/ Sort blocks /
6003	void
6004	ClpPackedMatrix3::sortBlocks(const ClpSimplex * model)
6005	{
6006	int * lookup = column_ + numberColumns_;
6007	for (int iBlock = `0`; iBlock < numberBlocks_; iBlock++) {
6008	blockStruct * block = block_ + iBlock;
6009	int numberInBlock = block->numberInBlock_;
6010	int nel = block->numberElements_;
6011	int * row = row_ + block->startElements_;
6012	double * element = element_ + block->startElements_;
6013	int * column = column_ + block->startIndices_;
6014	int lastPrice = `0`;
6015	int firstNotPrice = numberInBlock - `1`;
6016	while (lastPrice <= firstNotPrice) {
6017	// find first basic or fixed
6018	int iColumn = numberInBlock;
6019	for (; lastPrice <= firstNotPrice; lastPrice++) {
6020	iColumn = column[lastPrice];
6021	if (model->getColumnStatus(iColumn) == ClpSimplex::basic \|\|
6022	model->getColumnStatus(iColumn) == ClpSimplex::isFixed)
6023	break;
6024	}
6025	// find last non basic or fixed
6026	int jColumn = -`1`;
6027	for (; firstNotPrice > lastPrice; firstNotPrice--) {
6028	jColumn = column[firstNotPrice];
6029	if (model->getColumnStatus(jColumn) != ClpSimplex::basic &&
6030	model->getColumnStatus(jColumn) != ClpSimplex::isFixed)
6031	break;
6032	}
6033	if (firstNotPrice > lastPrice) {
6034	assert (column[lastPrice] == iColumn);
6035	assert (column[firstNotPrice] == jColumn);
6036	// need to swap
6037	column[firstNotPrice] = iColumn;
6038	lookup[iColumn] = firstNotPrice;
6039	column[lastPrice] = jColumn;
6040	lookup[jColumn] = lastPrice;
6041	double * elementA = element + lastPrice * nel;
6042	int * rowA = row + lastPrice * nel;
6043	double * elementB = element + firstNotPrice * nel;
6044	int * rowB = row + firstNotPrice * nel;
6045	for (int i = `0`; i < nel; i++) {
6046	int temp = rowA[i];
6047	double tempE = elementA[i];
6048	rowA[i] = rowB[i];
6049	elementA[i] = elementB[i];
6050	rowB[i] = temp;
6051	elementB[i] = tempE;
6052	}
6053	firstNotPrice--;
6054	lastPrice++;
6055	} else if (lastPrice == firstNotPrice) {
6056	// make sure correct side
6057	iColumn = column[lastPrice];
6058	if (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
6059	model->getColumnStatus(iColumn) != ClpSimplex::isFixed)
6060	lastPrice++;
6061	break;
6062	}
6063	}
6064	block->numberPrice_ = lastPrice;
6065	#ifndef NDEBUG
6066	// check
6067	int i;
6068	for (i = `0`; i < lastPrice; i++) {
6069	int iColumn = column[i];
6070	assert (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
6071	model->getColumnStatus(iColumn) != ClpSimplex::isFixed);
6072	assert (lookup[iColumn] == i);
6073	}
6074	for (; i < numberInBlock; i++) {
6075	int iColumn = column[i];
6076	assert (model->getColumnStatus(iColumn) == ClpSimplex::basic \|\|
6077	model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
6078	assert (lookup[iColumn] == i);
6079	}
6080	#endif
6081	}
6082	}
6083	// Swap one variable
6084	void
6085	ClpPackedMatrix3::swapOne(const ClpSimplex * model, const ClpPackedMatrix * matrix,
6086	int iColumn)
6087	{
6088	int * lookup = column_ + numberColumns_;
6089	// position in block
6090	int kA = lookup[iColumn];
6091	if (kA < `0`)
6092	return; // odd one
6093	// get matrix data pointers
6094	const CoinPackedMatrix * columnCopy = matrix->getPackedMatrix();
6095	//const int row = columnCopy->getIndices();*
6096	const CoinBigIndex * columnStart = columnCopy->getVectorStarts();
6097	const int * columnLength = columnCopy->getVectorLengths();
6098	const double * elementByColumn = columnCopy->getElements();
6099	CoinBigIndex start = columnStart[iColumn];
6100	int n = columnLength[iColumn];
6101	if (matrix->zeros()) {
6102	CoinBigIndex end = start + n;
6103	for (CoinBigIndex j = start; j < end; j++) {
6104	if(!elementByColumn[j])
6105	n--;
6106	}
6107	}
6108	// find block - could do binary search
6109	int iBlock = CoinMin(n, numberBlocks_) - `1`;
6110	while (block_[iBlock].numberElements_ != n)
6111	iBlock--;
6112	blockStruct * block = block_ + iBlock;
6113	int nel = block->numberElements_;
6114	int * row = row_ + block->startElements_;
6115	double * element = element_ + block->startElements_;
6116	int * column = column_ + block->startIndices_;
6117	assert (column[kA] == iColumn);
6118	bool moveUp = (model->getColumnStatus(iColumn) == ClpSimplex::basic \|\|
6119	model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
6120	int lastPrice = block->numberPrice_;
6121	int kB;
6122	if (moveUp) {
6123	// May already be in correct place (e.g. fixed basic leaving basis)
6124	if (kA >= lastPrice)
6125	return;
6126	kB = lastPrice - `1`;
6127	block->numberPrice_--;
6128	} else {
6129	assert (kA >= lastPrice);
6130	kB = lastPrice;
6131	block->numberPrice_++;
6132	}
6133	int jColumn = column[kB];
6134	column[kA] = jColumn;
6135	lookup[jColumn] = kA;
6136	column[kB] = iColumn;
6137	lookup[iColumn] = kB;
6138	double * elementA = element + kB * nel;
6139	int * rowA = row + kB * nel;
6140	double * elementB = element + kA * nel;
6141	int * rowB = row + kA * nel;
6142	int i;
6143	for (i = `0`; i < nel; i++) {
6144	int temp = rowA[i];
6145	double tempE = elementA[i];
6146	rowA[i] = rowB[i];
6147	elementA[i] = elementB[i];
6148	rowB[i] = temp;
6149	elementB[i] = tempE;
6150	}
6151	#ifndef NDEBUG
6152	// check
6153	for (i = `0`; i < block->numberPrice_; i++) {
6154	int iColumn = column[i];
6155	if (iColumn != model->sequenceIn() && iColumn != model->sequenceOut())
6156	assert (model->getColumnStatus(iColumn) != ClpSimplex::basic &&
6157	model->getColumnStatus(iColumn) != ClpSimplex::isFixed);
6158	assert (lookup[iColumn] == i);
6159	}
6160	int numberInBlock = block->numberInBlock_;
6161	for (; i < numberInBlock; i++) {
6162	int iColumn = column[i];
6163	if (iColumn != model->sequenceIn() && iColumn != model->sequenceOut())
6164	assert (model->getColumnStatus(iColumn) == ClpSimplex::basic \|\|
6165	model->getColumnStatus(iColumn) == ClpSimplex::isFixed);
6166	assert (lookup[iColumn] == i);
6167	}
6168	#endif
6169	}
6170	/ Return <code>x * -1 * A in <code>z</code>.*
6171	Note - x packed and z will be packed mode
6172	Squashes small elements and knows about ClpSimplex /*
6173	void
6174	ClpPackedMatrix3::transposeTimes(const ClpSimplex * model,
6175	const double * pi,
6176	CoinIndexedVector * output) const
6177	{
6178	int numberNonZero = `0`;
6179	int * index = output->getIndices();
6180	double * array = output->denseVector();
6181	double zeroTolerance = model->zeroTolerance();
6182	double value = `0.0`;
6183	CoinBigIndex j;
6184	int numberOdd = block_->startIndices_;
6185	if (numberOdd) {
6186	// A) as probably long may be worth unrolling
6187	CoinBigIndex end = start_[`1`];
6188	for (j = start_[`0`]; j < end; j++) {
6189	int iRow = row_[j];
6190	value += pi[iRow] * element_[j];
6191	}
6192	int iColumn;
6193	// int jColumn=column_[0];
6194
6195	for (iColumn = `0`; iColumn < numberOdd - `1`; iColumn++) {
6196	CoinBigIndex start = end;
6197	end = start_[iColumn+`2`];
6198	if (fabs(value) > zeroTolerance) {
6199	array[numberNonZero] = value;
6200	index[numberNonZero++] = column_[iColumn];
6201	//index[numberNonZero++]=jColumn;
6202	}
6203	// jColumn = column_[iColumn+1];
6204	value = `0.0`;
6205	//if (model->getColumnStatus(jColumn)!=ClpSimplex::basic) {
6206	for (j = start; j < end; j++) {
6207	int iRow = row_[j];
6208	value += pi[iRow] * element_[j];
6209	}
6210	//}
6211	}
6212	if (fabs(value) > zeroTolerance) {
6213	array[numberNonZero] = value;
6214	index[numberNonZero++] = column_[iColumn];
6215	//index[numberNonZero++]=jColumn;
6216	}
6217	}
6218	for (int iBlock = `0`; iBlock < numberBlocks_; iBlock++) {
6219	// B) Can sort so just do nonbasic (and nonfixed)
6220	// C) Can do two at a time (if so put odd one into start_)
6221	// D) can use switch
6222	blockStruct * block = block_ + iBlock;
6223	//int numberPrice = block->numberInBlock_;
6224	int numberPrice = block->numberPrice_;
6225	int nel = block->numberElements_;
6226	int * row = row_ + block->startElements_;
6227	double * element = element_ + block->startElements_;
6228	int * column = column_ + block->startIndices_;
6229	#if 0
6230	// two at a time
6231	if ((numberPrice & `1`) != `0`) {
6232	double value = `0.0`;
6233	int nel2 = nel;
6234	for (; nel2; nel2--) {
6235	int iRow = *row++;
6236	value += pi[iRow] * (*element++);
6237	}
6238	if (fabs(value) > zeroTolerance) {
6239	array[numberNonZero] = value;
6240	index[numberNonZero++] = *column;
6241	}
6242	column++;
6243	}
6244	numberPrice = numberPrice >> `1`;
6245	switch ((nel % `2`)) {
6246	// 2 k +0
6247	case `0`:
6248	for (; numberPrice; numberPrice--) {
6249	double value0 = `0.0`;
6250	double value1 = `0.0`;
6251	int nel2 = nel;
6252	for (; nel2; nel2--) {
6253	int iRow0 = *row;
6254	int iRow1 = *(row + nel);
6255	row++;
6256	double element0 = *element;
6257	double element1 = *(element + nel);
6258	element++;
6259	value0 += pi[iRow0] * element0;
6260	value1 += pi[iRow1] * element1;
6261	}
6262	row += nel;
6263	element += nel;
6264	if (fabs(value0) > zeroTolerance) {
6265	array[numberNonZero] = value0;
6266	index[numberNonZero++] = *column;
6267	}
6268	column++;
6269	if (fabs(value1) > zeroTolerance) {
6270	array[numberNonZero] = value1;
6271	index[numberNonZero++] = *column;
6272	}
6273	column++;
6274	}
6275	break;
6276	// 2 k +1
6277	case `1`:
6278	for (; numberPrice; numberPrice--) {
6279	double value0;
6280	double value1;
6281	int nel2 = nel - `1`;
6282	{
6283	int iRow0 = row[`0`];
6284	int iRow1 = row[nel];
6285	double pi0 = pi[iRow0];
6286	double pi1 = pi[iRow1];
6287	value0 = pi0 * element[`0`];
6288	value1 = pi1 * element[nel];
6289	row++;
6290	element++;
6291	}
6292	for (; nel2; nel2--) {
6293	int iRow0 = *row;
6294	int iRow1 = *(row + nel);
6295	row++;
6296	double element0 = *element;
6297	double element1 = *(element + nel);
6298	element++;
6299	value0 += pi[iRow0] * element0;
6300	value1 += pi[iRow1] * element1;
6301	}
6302	row += nel;
6303	element += nel;
6304	if (fabs(value0) > zeroTolerance) {
6305	array[numberNonZero] = value0;
6306	index[numberNonZero++] = *column;
6307	}
6308	column++;
6309	if (fabs(value1) > zeroTolerance) {
6310	array[numberNonZero] = value1;
6311	index[numberNonZero++] = *column;
6312	}
6313	column++;
6314	}
6315	break;
6316	}
6317	#else
6318	for (; numberPrice; numberPrice--) {
6319	double value = `0.0`;
6320	int nel2 = nel;
6321	for (; nel2; nel2--) {
6322	int iRow = *row++;
6323	value += pi[iRow] * (*element++);
6324	}
6325	if (fabs(value) > zeroTolerance) {
6326	array[numberNonZero] = value;
6327	index[numberNonZero++] = *column;
6328	}
6329	column++;
6330	}
6331	#endif
6332	}
6333	output->setNumElements(numberNonZero);
6334	}
6335	// Updates two arrays for steepest
6336	void
6337	ClpPackedMatrix3::transposeTimes2(const ClpSimplex * model,
6338	const double * pi, CoinIndexedVector * output,
6339	const double * piWeight,
6340	double referenceIn, double devex,
6341	// Array for exact devex to say what is in reference framework
6342	unsigned int * reference,
6343	double * weights, double scaleFactor)
6344	{
6345	int numberNonZero = `0`;
6346	int * index = output->getIndices();
6347	double * array = output->denseVector();
6348	double zeroTolerance = model->zeroTolerance();
6349	double value = `0.0`;
6350	bool killDjs = (scaleFactor == `0.0`);
6351	if (!scaleFactor)
6352	scaleFactor = `1.0`;
6353	int numberOdd = block_->startIndices_;
6354	int iColumn;
6355	CoinBigIndex end = start_[`0`];
6356	for (iColumn = `0`; iColumn < numberOdd; iColumn++) {
6357	CoinBigIndex start = end;
6358	CoinBigIndex j;
6359	int jColumn = column_[iColumn];
6360	end = start_[iColumn+`1`];
6361	value = `0.0`;
6362	if (model->getColumnStatus(jColumn) != ClpSimplex::basic) {
6363	for (j = start; j < end; j++) {
6364	int iRow = row_[j];
6365	value -= pi[iRow] * element_[j];
6366	}
6367	if (fabs(value) > zeroTolerance) {
6368	// and do other array
6369	double modification = `0.0`;
6370	for (j = start; j < end; j++) {
6371	int iRow = row_[j];
6372	modification += piWeight[iRow] * element_[j];
6373	}
6374	double thisWeight = weights[jColumn];
6375	double pivot = value * scaleFactor;
6376	double pivotSquared = pivot * pivot;
6377	thisWeight += pivotSquared * devex + pivot * modification;
6378	if (thisWeight < DEVEX_TRY_NORM) {
6379	if (referenceIn < `0.0`) {
6380	// steepest
6381	thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
6382	} else {
6383	// exact
6384	thisWeight = referenceIn * pivotSquared;
6385	if (reference(jColumn))
6386	thisWeight += `1.0`;
6387	thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
6388	}
6389	}
6390	weights[jColumn] = thisWeight;
6391	if (!killDjs) {
6392	array[numberNonZero] = value;
6393	index[numberNonZero++] = jColumn;
6394	}
6395	}
6396	}
6397	}
6398	for (int iBlock = `0`; iBlock < numberBlocks_; iBlock++) {
6399	// B) Can sort so just do nonbasic (and nonfixed)
6400	// C) Can do two at a time (if so put odd one into start_)
6401	// D) can use switch
6402	blockStruct * block = block_ + iBlock;
6403	//int numberPrice = block->numberInBlock_;
6404	int numberPrice = block->numberPrice_;
6405	int nel = block->numberElements_;
6406	int * row = row_ + block->startElements_;
6407	double * element = element_ + block->startElements_;
6408	int * column = column_ + block->startIndices_;
6409	for (; numberPrice; numberPrice--) {
6410	double value = `0.0`;
6411	int nel2 = nel;
6412	for (; nel2; nel2--) {
6413	int iRow = *row++;
6414	value -= pi[iRow] * (*element++);
6415	}
6416	if (fabs(value) > zeroTolerance) {
6417	int jColumn = *column;
6418	// back to beginning
6419	row -= nel;
6420	element -= nel;
6421	// and do other array
6422	double modification = `0.0`;
6423	nel2 = nel;
6424	for (; nel2; nel2--) {
6425	int iRow = *row++;
6426	modification += piWeight[iRow] * (*element++);
6427	}
6428	double thisWeight = weights[jColumn];
6429	double pivot = value * scaleFactor;
6430	double pivotSquared = pivot * pivot;
6431	thisWeight += pivotSquared * devex + pivot * modification;
6432	if (thisWeight < DEVEX_TRY_NORM) {
6433	if (referenceIn < `0.0`) {
6434	// steepest
6435	thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
6436	} else {
6437	// exact
6438	thisWeight = referenceIn * pivotSquared;
6439	if (reference(jColumn))
6440	thisWeight += `1.0`;
6441	thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
6442	}
6443	}
6444	weights[jColumn] = thisWeight;
6445	if (!killDjs) {
6446	array[numberNonZero] = value;
6447	index[numberNonZero++] = jColumn;
6448	}
6449	}
6450	column++;
6451	}
6452	}
6453	output->setNumElements(numberNonZero);
6454	output->setPackedMode(true);
6455	}
6456	#if COIN_LONG_WORK
6457	// For long double versions
6458	void
6459	ClpPackedMatrix::times(CoinWorkDouble scalar,
6460	const CoinWorkDouble * x, CoinWorkDouble * y) const
6461	{
6462	int iRow, iColumn;
6463	// get matrix data pointers
6464	const int * row = matrix_->getIndices();
6465	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
6466	const double * elementByColumn = matrix_->getElements();
6467	//memset(y,0,matrix_->getNumRows()sizeof(double));*
6468	assert (((flags_ & `2`) != `0`) == matrix_->hasGaps());
6469	if (!(flags_ & `2`)) {
6470	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
6471	CoinBigIndex j;
6472	CoinWorkDouble value = x[iColumn];
6473	if (value) {
6474	CoinBigIndex start = columnStart[iColumn];
6475	CoinBigIndex end = columnStart[iColumn+`1`];
6476	value *= scalar;
6477	for (j = start; j < end; j++) {
6478	iRow = row[j];
6479	y[iRow] += value * elementByColumn[j];
6480	}
6481	}
6482	}
6483	} else {
6484	const int * columnLength = matrix_->getVectorLengths();
6485	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
6486	CoinBigIndex j;
6487	CoinWorkDouble value = x[iColumn];
6488	if (value) {
6489	CoinBigIndex start = columnStart[iColumn];
6490	CoinBigIndex end = start + columnLength[iColumn];
6491	value *= scalar;
6492	for (j = start; j < end; j++) {
6493	iRow = row[j];
6494	y[iRow] += value * elementByColumn[j];
6495	}
6496	}
6497	}
6498	}
6499	}
6500	void
6501	ClpPackedMatrix::transposeTimes(CoinWorkDouble scalar,
6502	const CoinWorkDouble * x, CoinWorkDouble * y) const
6503	{
6504	int iColumn;
6505	// get matrix data pointers
6506	const int * row = matrix_->getIndices();
6507	const CoinBigIndex * columnStart = matrix_->getVectorStarts();
6508	const double * elementByColumn = matrix_->getElements();
6509	if (!(flags_ & `2`)) {
6510	if (scalar == -`1.0`) {
6511	CoinBigIndex start = columnStart[`0`];
6512	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
6513	CoinBigIndex j;
6514	CoinBigIndex next = columnStart[iColumn+`1`];
6515	CoinWorkDouble value = y[iColumn];
6516	for (j = start; j < next; j++) {
6517	int jRow = row[j];
6518	value -= x[jRow] * elementByColumn[j];
6519	}
6520	start = next;
6521	y[iColumn] = value;
6522	}
6523	} else {
6524	CoinBigIndex start = columnStart[`0`];
6525	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
6526	CoinBigIndex j;
6527	CoinBigIndex next = columnStart[iColumn+`1`];
6528	CoinWorkDouble value = `0.0`;
6529	for (j = start; j < next; j++) {
6530	int jRow = row[j];
6531	value += x[jRow] * elementByColumn[j];
6532	}
6533	start = next;
6534	y[iColumn] += value * scalar;
6535	}
6536	}
6537	} else {
6538	const int * columnLength = matrix_->getVectorLengths();
6539	for (iColumn = `0`; iColumn < numberActiveColumns_; iColumn++) {
6540	CoinBigIndex j;
6541	CoinWorkDouble value = `0.0`;
6542	CoinBigIndex start = columnStart[iColumn];
6543	CoinBigIndex end = start + columnLength[iColumn];
6544	for (j = start; j < end; j++) {
6545	int jRow = row[j];
6546	value += x[jRow] * elementByColumn[j];
6547	}
6548	y[iColumn] += value * scalar;
6549	}
6550	}
6551	}
6552	#endif
6553	#ifdef CLP_ALL_ONE_FILE
6554	#undef reference
6555	#endif
6556

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