CoinPackedMatrix.cpp source code [OGDF/src/coin/CoinUtils/CoinPackedMatrix.cpp]

1	/ $Id: CoinPackedMatrix.cpp 1539 2012-06-28 10:26:15Z forrest $ /
2	// Copyright (C) 2000, International Business Machines
3	// Corporation and others. All Rights Reserved.
4	// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6	#include "CoinUtilsConfig.h"
7
8	#include <algorithm>
9	#include <numeric>
10	#include <cassert>
11	#include <cstdio>
12	#include <cmath>
13	#include <iostream>
14
15	#include "CoinPragma.hpp"
16	#include "CoinSort.hpp"
17	#include "CoinHelperFunctions.hpp"
18	#ifndef CLP_NO_VECTOR
19	#include "CoinPackedVectorBase.hpp"
20	#endif
21	#include "CoinFloatEqual.hpp"
22	#include "CoinPackedMatrix.hpp"
23
24	#if !defined(COIN_COINUTILS_CHECKLEVEL)
25	#define COIN_COINUTILS_CHECKLEVEL 0
26	#endif
27
28	//#############################################################################
29	// T must be an integral type (int, CoinBigIndex, etc.)
30	template <typename T>
31	static inline T
32	CoinLengthWithExtra(T len, double extraGap)
33	{
34	return static_cast<T>(ceil(len * (`1` + extraGap)));
35	}
36
37	//#############################################################################
38
39	static inline void
40	CoinTestSortedIndexSet(const int num, const int * sorted, const int maxEntry,
41	const char * testingMethod)
42	{
43	if (sorted[`0`] < `0` \|\| sorted[num-`1`] >= maxEntry)
44	throw CoinError ("bad index", testingMethod, "CoinPackedMatrix");
45	if (std::adjacent_find(sorted, sorted + num) != sorted + num)
46	throw CoinError ("duplicate index", testingMethod, "CoinPackedMatrix");
47	}
48
49	//-----------------------------------------------------------------------------
50
51	static inline int *
52	CoinTestIndexSet(const int numDel, const int * indDel, const int maxEntry,
53	const char * testingMethod)
54	{
55	if (! CoinIsSorted(indDel, indDel + numDel)) {
56	// if not sorted then sort it, test for consistency and return a pointer
57	// to the sorted array
58	int * sorted = new int[numDel];
59	CoinMemcpyN(indDel, numDel, sorted);
60	std::sort(sorted, sorted + numDel);
61	CoinTestSortedIndexSet(numDel, sorted, maxEntry, testingMethod);
62	return sorted;
63	}
64
65	// Otherwise it's already sorted, so just test for consistency and return a
66	// 0 pointer.
67	CoinTestSortedIndexSet(numDel, indDel, maxEntry, testingMethod);
68	return `0`;
69	}
70
71	//#############################################################################
72
73	void
74	CoinPackedMatrix::reserve(const int newMaxMajorDim, const CoinBigIndex newMaxSize,
75	bool create)
76	{
77	if (newMaxMajorDim > maxMajorDim_) {
78	maxMajorDim_ = newMaxMajorDim;
79	int * oldlength = length_;
80	CoinBigIndex * oldstart = start_;
81	length_ = new int[newMaxMajorDim];
82	start_ = new CoinBigIndex[newMaxMajorDim+`1`];
83	start_[`0`]=`0`;
84	if (majorDim_ > `0`) {
85	CoinMemcpyN(oldlength, majorDim_, length_);
86	CoinMemcpyN(oldstart, majorDim_ + `1`, start_);
87	}
88	if (create) {
89	// create empty vectors
90	CoinFillN(length_+majorDim_,maxMajorDim_-majorDim_,`0`);
91	CoinFillN(start_+majorDim_+`1`,maxMajorDim_-majorDim_,`0`);
92	majorDim_=maxMajorDim_;
93	}
94	delete[] oldlength;
95	delete[] oldstart;
96	}
97	if (newMaxSize > maxSize_) {
98	maxSize_ = newMaxSize;
99	int * oldind = index_;
100	double * oldelem = element_;
101	index_ = new int[newMaxSize];
102	element_ = new double[newMaxSize];
103	for (int i = majorDim_ - `1`; i >= `0`; --i) {
104	CoinMemcpyN(oldind+start_[i], length_[i], index_+start_[i]);
105	CoinMemcpyN(oldelem+start_[i], length_[i], element_+start_[i]);
106	}
107	delete[] oldind;
108	delete[] oldelem;
109	}
110	}
111
112	//-----------------------------------------------------------------------------
113
114	void
115	CoinPackedMatrix::clear()
116	{
117	majorDim_ = `0`;
118	minorDim_ = `0`;
119	size_ = `0`;
120	}
121
122	//#############################################################################
123	//#############################################################################
124
125	void
126	CoinPackedMatrix::setDimensions(int newnumrows, int newnumcols)
127	{
128	const int numrows = getNumRows();
129	if (newnumrows < `0`)
130	newnumrows = numrows;
131	if (newnumrows < numrows)
132	throw CoinError ("Bad new rownum (less than current)",
133	"setDimensions", "CoinPackedMatrix");
134
135	const int numcols = getNumCols();
136	if (newnumcols < `0`)
137	newnumcols = numcols;
138	if (newnumcols < numcols)
139	throw CoinError ("Bad new colnum (less than current)",
140	"setDimensions", "CoinPackedMatrix");
141
142	int numplus = `0`;
143	if (isColOrdered()) {
144	minorDim_ = newnumrows;
145	numplus = newnumcols - numcols;
146	} else {
147	minorDim_ = newnumcols;
148	numplus = newnumrows - numrows;
149	}
150	if (numplus > `0`) {
151	int* lengths = new int[numplus];
152	CoinZeroN(lengths, numplus);
153	resizeForAddingMajorVectors(numplus, lengths);
154	delete[] lengths;
155	majorDim_ += numplus; //forgot to change majorDim_
156	}
157
158	}
159
160	//-----------------------------------------------------------------------------
161
162	void
163	CoinPackedMatrix::setExtraGap(const double newGap)
164	{
165	if (newGap < `0`)
166	throw CoinError ("negative new extra gap",
167	"setExtraGap", "CoinPackedMatrix");
168	extraGap_ = newGap;
169	}
170
171	//-----------------------------------------------------------------------------
172
173	void
174	CoinPackedMatrix::setExtraMajor(const double newMajor)
175	{
176	if (newMajor < `0`)
177	throw CoinError ("negative new extra major",
178	"setExtraMajor", "CoinPackedMatrix");
179	extraMajor_ = newMajor;
180	}
181
182	//#############################################################################
183	#ifndef CLP_NO_VECTOR
184	void
185	CoinPackedMatrix::appendCol(const CoinPackedVectorBase& vec)
186	{
187	if (colOrdered_)
188	appendMajorVector(vec);
189	else
190	appendMinorVector(vec);
191	}
192	#endif
193	//-----------------------------------------------------------------------------
194
195	void
196	CoinPackedMatrix::appendCol(const int vecsize,
197	const int *vecind,
198	const double *vecelem)
199	{
200	if (colOrdered_)
201	appendMajorVector(vecsize, vecind, vecelem);
202	else
203	appendMinorVector(vecsize, vecind, vecelem);
204	}
205
206	//-----------------------------------------------------------------------------
207	#ifndef CLP_NO_VECTOR
208	void
209	CoinPackedMatrix::appendCols(const int numcols,
210	const CoinPackedVectorBase * const * cols)
211	{
212	if (colOrdered_)
213	appendMajorVectors(numcols, cols);
214	else
215	appendMinorVectors(numcols, cols);
216	}
217	#endif
218	//-----------------------------------------------------------------------------
219
220	int
221	CoinPackedMatrix::appendCols(const int numcols,
222	const CoinBigIndex * columnStarts, const int * row,
223	const double * element, int numberRows)
224	{
225	int numberErrors;
226	if (colOrdered_) {
227	numberErrors=appendMajor(numcols, columnStarts, row, element, numberRows);
228	} else {
229	numberErrors=appendMinor(numcols, columnStarts, row, element, numberRows);
230	}
231	return numberErrors;
232	}
233	//-----------------------------------------------------------------------------
234	#ifndef CLP_NO_VECTOR
235	void
236	CoinPackedMatrix::appendRow(const CoinPackedVectorBase& vec)
237	{
238	if (colOrdered_)
239	appendMinorVector(vec);
240	else
241	appendMajorVector(vec);
242	}
243	#endif
244	//-----------------------------------------------------------------------------
245
246	void
247	CoinPackedMatrix::appendRow(const int vecsize,
248	const int *vecind,
249	const double *vecelem)
250	{
251	if (colOrdered_)
252	appendMinorVector(vecsize, vecind, vecelem);
253	else
254	appendMajorVector(vecsize, vecind, vecelem);
255	}
256
257	//-----------------------------------------------------------------------------
258	#ifndef CLP_NO_VECTOR
259	void
260	CoinPackedMatrix::appendRows(const int numrows,
261	const CoinPackedVectorBase * const * rows)
262	{
263	if (colOrdered_) {
264	// make sure enough columns
265	if (numrows == `0`)
266	return;
267
268	int i;
269	int maxDim=-`1`;
270	for (i = numrows - `1`; i >= `0`; --i) {
271	const int vecsize = rows[i]->getNumElements();
272	const int* vecind = rows[i]->getIndices();
273	for (int j = vecsize - `1`; j >= `0`; --j)
274	maxDim = CoinMax(maxDim,vecind[j]);
275	}
276	maxDim++;
277	if (maxDim>majorDim_) {
278	setDimensions(minorDim_,maxDim);
279	//int nAdd=maxDim-majorDim_;
280	//int length = new int[nAdd];*
281	//memset(length,0,nAddsizeof(int));*
282	//resizeForAddingMajorVectors(nAdd,length);
283	//delete [] length;
284	}
285	appendMinorVectors(numrows, rows);
286	} else {
287	appendMajorVectors(numrows, rows);
288	}
289	}
290	#endif
291	//-----------------------------------------------------------------------------
292
293	int
294	CoinPackedMatrix::appendRows(const int numrows,
295	const CoinBigIndex * rowStarts, const int * column,
296	const double * element, int numberColumns)
297	{
298	int numberErrors;
299	if (colOrdered_) {
300	numberErrors=appendMinor(numrows, rowStarts, column, element, numberColumns);
301	} else {
302	numberErrors=appendMajor(numrows, rowStarts, column, element, numberColumns);
303	}
304	return numberErrors;
305	}
306
307	//#############################################################################
308
309	void
310	CoinPackedMatrix::rightAppendPackedMatrix(const CoinPackedMatrix& matrix)
311	{
312	if (colOrdered_) {
313	if (matrix.colOrdered_) {
314	majorAppendSameOrdered(matrix);
315	} else {
316	majorAppendOrthoOrdered(matrix);
317	}
318	} else {
319	if (matrix.colOrdered_) {
320	minorAppendOrthoOrdered(matrix);
321	} else {
322	minorAppendSameOrdered(matrix);
323	}
324	}
325	}
326
327	//-----------------------------------------------------------------------------
328
329	void
330	CoinPackedMatrix::bottomAppendPackedMatrix(const CoinPackedMatrix& matrix)
331	{
332	if (colOrdered_) {
333	if (matrix.colOrdered_) {
334	minorAppendSameOrdered(matrix);
335	} else {
336	minorAppendOrthoOrdered(matrix);
337	}
338	} else {
339	if (matrix.colOrdered_) {
340	majorAppendOrthoOrdered(matrix);
341	} else {
342	majorAppendSameOrdered(matrix);
343	}
344	}
345	}
346
347	//#############################################################################
348
349	void
350	CoinPackedMatrix::deleteCols(const int numDel, const int * indDel)
351	{
352	if (numDel) {
353	if (colOrdered_)
354	deleteMajorVectors(numDel, indDel);
355	else
356	deleteMinorVectors(numDel, indDel);
357	}
358	}
359
360	//-----------------------------------------------------------------------------
361
362	void
363	CoinPackedMatrix::deleteRows(const int numDel, const int * indDel)
364	{
365	if (numDel) {
366	if (colOrdered_)
367	deleteMinorVectors(numDel, indDel);
368	else
369	deleteMajorVectors(numDel, indDel);
370	}
371	}
372
373	//#############################################################################
374	/ Replace the elements of a vector. The indices remain the same.*
375	At most the number specified will be replaced.
376	The index is between 0 and major dimension of matrix /*
377	void
378	CoinPackedMatrix::replaceVector(const int index,
379	const int numReplace,
380	const double * newElements)
381	{
382	if (index >= `0` && index < majorDim_) {
383	int length = (length_[index] < numReplace) ? length_[index] : numReplace;
384	CoinMemcpyN(newElements, length, element_ + start_[index]);
385	} else {
386	#ifdef COIN_DEBUG
387	throw CoinError("bad index", "replaceVector", "CoinPackedMatrix");
388	#endif
389	}
390	}
391	/ Modify one element of packed matrix. An element may be added.*
392	If the new element is zero it will be deleted unless
393	keepZero true /*
394	void
395	CoinPackedMatrix::modifyCoefficient(int row, int column, double newElement,
396	bool keepZero)
397	{
398	int minorIndex,majorIndex;
399	if (colOrdered_) {
400	majorIndex=column;
401	minorIndex=row;
402	} else {
403	minorIndex=column;
404	majorIndex=row;
405	}
406	if (majorIndex >= `0` && majorIndex < majorDim_) {
407	if (minorIndex >= `0` && minorIndex < minorDim_) {
408	CoinBigIndex j;
409	CoinBigIndex end=start_[majorIndex]+length_[majorIndex];
410	for (j=start_[majorIndex];j<end;j++) {
411	if (minorIndex==index_[j]) {
412	// replacement
413	if (newElement\|\|keepZero) {
414	element_[j]=newElement;
415	} else {
416	// pack down and return
417	length_[majorIndex]--;
418	end--;
419	size_--;
420	for (;j<end;j++) {
421	element_[j]=element_[j+`1`];
422	index_[j]=index_[j+`1`];
423	}
424	}
425	return;
426	}
427	}
428	if (j==end&&(newElement\|\|keepZero)) {
429	// we need to insert - keep in minor order if possible
430	if (end>=start_[majorIndex+`1`]) {
431	int * addedEntries = new int[majorDim_];
432	memset(addedEntries, `0`, majorDim_ * sizeof(int));
433	addedEntries[majorIndex] = `1`;
434	resizeForAddingMinorVectors(addedEntries);
435	delete[] addedEntries;
436	}
437	// So where to insert? We're just going to assume that the entries
438	// in the major vector are in increasing order, so we'll insert the
439	// new entry to the last place we can
440	const CoinBigIndex start = start_[majorIndex];
441	end = start_[majorIndex]+length_[majorIndex]; // recalculate end
442	for (j = end - `1`; j >= start; --j) {
443	if (index_[j] < minorIndex)
444	break;
445	index_[j+`1`] = index_[j];
446	element_[j+`1`] = element_[j];
447	}
448	++j;
449	index_[j] = minorIndex;
450	element_[j] = newElement;
451	size_++;
452	length_[majorIndex]++;
453	}
454	} else {
455	#ifdef COIN_DEBUG
456	throw CoinError("bad minor index", "modifyCoefficient",
457	"CoinPackedMatrix");
458	#endif
459	}
460	} else {
461	#ifdef COIN_DEBUG
462	throw CoinError("bad major index", "modifyCoefficient",
463	"CoinPackedMatrix");
464	#endif
465	}
466	}
467	/ Return one element of packed matrix.*
468	This works for either ordering
469	If it is not present will return 0.0 /*
470	double
471	CoinPackedMatrix::getCoefficient(int row, int column) const
472	{
473	int minorIndex,majorIndex;
474	if (colOrdered_) {
475	majorIndex=column;
476	minorIndex=row;
477	} else {
478	minorIndex=column;
479	majorIndex=row;
480	}
481	double value=`0.0`;
482	if (majorIndex >= `0` && majorIndex < majorDim_) {
483	if (minorIndex >= `0` && minorIndex < minorDim_) {
484	CoinBigIndex j;
485	CoinBigIndex end=start_[majorIndex]+length_[majorIndex];
486	for (j=start_[majorIndex];j<end;j++) {
487	if (minorIndex==index_[j]) {
488	value = element_[j];
489	break;
490	}
491	}
492	} else {
493	#ifdef COIN_DEBUG
494	throw CoinError("bad minor index", "modifyCoefficient",
495	"CoinPackedMatrix");
496	#endif
497	}
498	} else {
499	#ifdef COIN_DEBUG
500	throw CoinError("bad major index", "modifyCoefficient",
501	"CoinPackedMatrix");
502	#endif
503	}
504	return value;
505	}
506
507	//#############################################################################
508	/ Eliminate all elements in matrix whose*
509	absolute value is less than threshold.
510	The column starts are not affected. Returns number of elements
511	eliminated. Elements eliminated are at end of each vector
512	*/
513	int
514	CoinPackedMatrix::compress(double threshold)
515	{
516	CoinBigIndex numberEliminated =`0`;
517	// space for eliminated
518	int * eliminatedIndex = new int[minorDim_];
519	double * eliminatedElement = new double[minorDim_];
520	int i;
521	for (i=`0`;i<majorDim_;i++) {
522	int length = length_[i];
523	CoinBigIndex k=start_[i];
524	int kbad=`0`;
525	CoinBigIndex j;
526	for (j=start_[i];j<start_[i]+length;j++) {
527	if (fabs(element_[j])>=threshold) {
528	element_[k]=element_[j];
529	index_[k++]=index_[j];
530	} else {
531	eliminatedElement[kbad]=element_[j];
532	eliminatedIndex[kbad++]=index_[j];
533	}
534	}
535	if (kbad) {
536	numberEliminated += kbad;
537	length_[i] = k-start_[i];
538	memcpy(index_+k,eliminatedIndex,kbad*sizeof(int));
539	memcpy(element_+k,eliminatedElement,kbad*sizeof(double));
540	}
541	}
542	size_ -= numberEliminated;
543	delete [] eliminatedIndex;
544	delete [] eliminatedElement;
545	return numberEliminated;
546	}
547	//#############################################################################
548	/ Eliminate all elements in matrix whose*
549	absolute value is less than threshold.ALSO removes duplicates
550	The column starts are not affected. Returns number of elements
551	eliminated.
552	*/
553	int
554	CoinPackedMatrix::eliminateDuplicates(double threshold)
555	{
556	CoinBigIndex numberEliminated =`0`;
557	// space for eliminated
558	int * mark = new int [minorDim_];
559	int i;
560	for (i=`0`;i<minorDim_;i++)
561	mark[i]=-`1`;
562	for (i=`0`;i<majorDim_;i++) {
563	CoinBigIndex k=start_[i];
564	CoinBigIndex end = k+length_[i];
565	CoinBigIndex j;
566	for (j=k;j<end;j++) {
567	int index = index_[j];
568	if (mark[index]==-`1`) {
569	mark[index]=j;
570	} else {
571	// duplicate
572	int jj = mark[index];
573	element_[jj] += element_[j];
574	element_[j]=`0.0`;
575	}
576	}
577	for (j=k;j<end;j++) {
578	int index = index_[j];
579	mark[index]=-`1`;
580	if (fabs(element_[j])>=threshold) {
581	element_[k]=element_[j];
582	index_[k++]=index_[j];
583	}
584	}
585	numberEliminated += end-k;
586	length_[i] = k-start_[i];
587	}
588	size_ -= numberEliminated;
589	delete [] mark;
590	return numberEliminated;
591	}
592	//#############################################################################
593
594	void
595	CoinPackedMatrix::removeGaps(double removeValue)
596	{
597	if (removeValue<`0.0`) {
598	if (size_<start_[majorDim_]) {
599	#if 1
600	// Small copies so faster to do simply
601	int i;
602	CoinBigIndex size=`0`;
603	for (i = `1`; i < majorDim_+`1`; ++i) {
604	const CoinBigIndex si = start_[i];
605	size += length_[i-`1`];
606	if (si>size)
607	break;
608	}
609	for (; i < majorDim_; ++i) {
610	const CoinBigIndex si = start_[i];
611	const int li = length_[i];
612	start_[i] = size;
613	for (CoinBigIndex j=si;j<si+li;j++) {
614	assert (size<size_);
615	index_[size]=index_[j];
616	element_[size++]=element_[j];
617	}
618	}
619	assert (size==size_);
620	start_[majorDim_] = size;
621	for (i=`0`; i < majorDim_; ++i) {
622	assert (start_[i+`1`]==start_[i]+length_[i]);
623	}
624	#else
625	for (int i = `1`; i < majorDim_; ++i) {
626	const CoinBigIndex si = start_[i];
627	const int li = length_[i];
628	start_[i] = start_[i-`1`] + length_[i-`1`];
629	CoinCopy(index_ + si, index_ + (si + li), index_ + start_[i]);
630	CoinCopy(element_ + si, element_ + (si + li), element_ + start_[i]);
631
632	}
633	start_[majorDim_] = size_;
634	#endif
635	} else {
636	#ifndef NDEBUG
637	for (int i = `1`; i < majorDim_; ++i) {
638	assert (start_[i] == start_[i-`1`] + length_[i-`1`]);
639	}
640	assert(start_[majorDim_] == size_);
641	#endif
642	}
643	} else {
644	CoinBigIndex put=`0`;
645	assert (!start_[`0`]);
646	CoinBigIndex start = `0`;
647	for (int i = `0`; i < majorDim_; ++i) {
648	const CoinBigIndex si = start;
649	start = start_[i+`1`];
650	const int li = length_[i];
651	for (CoinBigIndex j = si;j<si+li;j++) {
652	double value = element_[j];
653	if (fabs(value)>removeValue) {
654	index_[put]=index_[j];
655	element_[put++]=value;
656	}
657	}
658	length_[i]=put-start_[i];
659	start_[i+`1`] = put;
660	}
661	size_ = put;
662	}
663	}
664
665	//#############################################################################
666
667	/ Really clean up matrix.*
668	a) eliminate all duplicate AND small elements in matrix
669	b) remove all gaps and set extraGap_ and extraMajor_ to 0.0
670	c) reallocate arrays and make max lengths equal to lengths
671	d) orders elements
672	returns number of elements eliminated
673	*/
674	int
675	CoinPackedMatrix::cleanMatrix(double threshold)
676	{
677	if (!majorDim_) {
678	extraGap_=`0.0`;
679	extraMajor_=`0.0`;
680	return `0`;
681	}
682	CoinBigIndex numberEliminated =`0`;
683	// space for eliminated
684	int * mark = new int [minorDim_];
685	int i;
686	for (i=`0`;i<minorDim_;i++)
687	mark[i]=-`1`;
688	CoinBigIndex n = `0`;
689	for (i=`0`;i<majorDim_;i++) {
690	CoinBigIndex k=start_[i];
691	start_[i]=n;
692	CoinBigIndex end = k+length_[i];
693	CoinBigIndex j;
694	for (j=k;j<end;j++) {
695	int index = index_[j];
696	if (mark[index]==-`1`) {
697	mark[index]=j;
698	} else {
699	// duplicate
700	int jj = mark[index];
701	element_[jj] += element_[j];
702	element_[j]=`0.0`;
703	}
704	}
705	for (j=k;j<end;j++) {
706	int index = index_[j];
707	mark[index]=-`1`;
708	if (fabs(element_[j])>=threshold) {
709	element_[n]=element_[j];
710	index_[n++]=index_[j];
711	k++;
712	}
713	}
714	numberEliminated += end-k;
715	length_[i] = n-start_[i];
716	// sort
717	CoinSort_2(index_+start_[i],index_+n,element_+start_[i]);
718	}
719	start_[majorDim_]=n;
720	size_ -= numberEliminated;
721	assert (n==size_);
722	delete [] mark;
723	extraGap_=`0.0`;
724	extraMajor_=`0.0`;
725	maxMajorDim_=majorDim_;
726	maxSize_=size_;
727	// Now reallocate - do smallest ones first
728	int * temp = CoinCopyOfArray(length_,majorDim_);
729	delete [] length_;
730	length_ = temp;
731	CoinBigIndex * temp2 = CoinCopyOfArray(start_,majorDim_+`1`);
732	delete [] start_;
733	start_ = temp2;
734	temp = CoinCopyOfArray(index_,size_);
735	delete [] index_;
736	index_ = temp;
737	double * temp3 = CoinCopyOfArray(element_,size_);
738	delete [] element_;
739	element_ = temp3;
740	return numberEliminated;
741	}
742
743	//#############################################################################
744
745	void
746	CoinPackedMatrix::submatrixOf(const CoinPackedMatrix& matrix,
747	const int numMajor, const int * indMajor)
748	{
749	int i;
750	int* sortedIndPtr = CoinTestIndexSet(numMajor, indMajor, matrix.majorDim_,
751	"submatrixOf");
752	const int * sortedInd = sortedIndPtr == `0` ? indMajor : sortedIndPtr;
753
754	gutsOfDestructor();
755
756	// Count how many nonzeros there'll be
757	CoinBigIndex nzcnt = `0`;
758	const int* length = matrix.getVectorLengths();
759	for (i = `0`; i < numMajor; ++i) {
760	nzcnt += length[sortedInd[i]];
761	}
762
763	colOrdered_ = matrix.colOrdered_;
764	maxMajorDim_ = int(numMajor * (`1`+extraMajor_) + `1`);
765	maxSize_ = static_cast<CoinBigIndex> (nzcnt * (`1`+extraMajor_) * (`1`+extraGap_) + `100`);
766	length_ = new int[maxMajorDim_];
767	start_ = new CoinBigIndex[maxMajorDim_+`1`];
768	start_[`0`]=`0`;
769	index_ = new int[maxSize_];
770	element_ = new double[maxSize_];
771	majorDim_ = `0`;
772	minorDim_ = matrix.minorDim_;
773	size_ = `0`;
774	#ifdef CLP_NO_VECTOR
775	for (i = `0`; i < numMajor; ++i) {
776	int j = sortedInd[i];
777	CoinBigIndex start = matrix.start_[j];
778	appendMajorVector(matrix.length_[j],matrix.index_+start,matrix.element_+start);
779	}
780	#else
781	for (i = `0`; i < numMajor; ++i) {
782	const CoinShallowPackedVector reqdBySunCC = matrix.getVector(sortedInd[i]) ;
783	appendMajorVector(reqdBySunCC);
784	}
785	#endif
786
787	delete[] sortedIndPtr;
788	}
789
790	//#############################################################################
791
792	void
793	CoinPackedMatrix::submatrixOfWithDuplicates(const CoinPackedMatrix& matrix,
794	const int numMajor, const int * indMajor)
795	{
796	int i;
797	// we allow duplicates - can be useful
798	#ifndef NDEBUG
799	for (i=`0`; i<numMajor;i++) {
800	if (indMajor[i]<`0`\|\|indMajor[i]>=matrix.majorDim_)
801	throw CoinError("bad index", "submatrixOfWithDuplicates", "CoinPackedMatrix");
802	}
803	#endif
804	gutsOfDestructor();
805	// Get rid of gaps
806	extraMajor_ = `0`;
807	extraGap_ = `0`;
808	colOrdered_ = matrix.colOrdered_;
809	maxMajorDim_ = numMajor ;
810
811	const int* length = matrix.getVectorLengths();
812	length_ = new int[maxMajorDim_];
813	start_ = new CoinBigIndex[maxMajorDim_+`1`];
814	// Count how many nonzeros there'll be
815	CoinBigIndex nzcnt = `0`;
816	for (i = `0`; i < maxMajorDim_; ++i) {
817	start_[i]=nzcnt;
818	int thisLength = length[indMajor[i]];
819	nzcnt += thisLength;
820	length_[i]=thisLength;
821	}
822	start_[maxMajorDim_]=nzcnt;
823	maxSize_ = nzcnt ;
824	index_ = new int[maxSize_];
825	element_ = new double[maxSize_];
826	majorDim_ = maxMajorDim_;
827	minorDim_ = matrix.minorDim_;
828	size_ = `0`;
829	const CoinBigIndex * startOld = matrix.start_;
830	const double * elementOld = matrix.element_;
831	const int * indexOld = matrix.index_;
832	for (i = `0`; i < maxMajorDim_; ++i) {
833	int j = indMajor[i];
834	CoinBigIndex start = startOld[j];
835	int thisLength = length_[i];
836	const double * element = elementOld+start;
837	const int * index = indexOld+start;
838	for (int j=`0`;j<thisLength;j++) {
839	element_[size_] = element[j];
840	index_[size_++] = index[j];
841	}
842	}
843	}
844
845	//#############################################################################
846
847	void
848	CoinPackedMatrix::copyOf(const CoinPackedMatrix& rhs)
849	{
850	if (this != &rhs) {
851	gutsOfDestructor();
852	gutsOfCopyOf(rhs.colOrdered_,
853	rhs.minorDim_, rhs.majorDim_, rhs.size_,
854	rhs.element_, rhs.index_, rhs.start_, rhs.length_,
855	rhs.extraMajor_, rhs.extraGap_);
856	}
857	}
858
859	//-----------------------------------------------------------------------------
860
861	void
862	CoinPackedMatrix::copyOf(const bool colordered,
863	const int minor, const int major,
864	const CoinBigIndex numels,
865	const double * elem, const int * ind,
866	const CoinBigIndex * start, const int * len,
867	const double extraMajor, const double extraGap)
868	{
869	gutsOfDestructor();
870	gutsOfCopyOf(colordered, minor, major, numels, elem, ind, start, len,
871	extraMajor, extraGap);
872	}
873	//#############################################################################
874	/ Copy method. This method makes an exact replica of the argument,*
875	including the extra space parameters.
876	If there is room it will re-use arrays /*
877	void
878	CoinPackedMatrix::copyReuseArrays(const CoinPackedMatrix& rhs)
879	{
880	assert (colOrdered_==rhs.colOrdered_);
881	if (maxMajorDim_>=rhs.majorDim_&&maxSize_>=rhs.size_) {
882	majorDim_ = rhs.majorDim_;
883	minorDim_ = rhs.minorDim_;
884	size_ = rhs.size_;
885	extraGap_ = rhs.extraGap_;
886	extraMajor_ = rhs.extraMajor_;
887	CoinMemcpyN(rhs.length_, majorDim_,length_);
888	CoinMemcpyN(rhs.start_, majorDim_+`1`,start_);
889	if (size_==start_[majorDim_]) {
890	CoinMemcpyN(rhs.index_ , size_, index_);
891	CoinMemcpyN(rhs.element_ , size_, element_);
892	} else {
893	// we can't just simply memcpy these content over, because that can
894	// upset memory debuggers like purify if there were gaps and those gaps
895	// were uninitialized memory blocks
896	for (int i = majorDim_ - `1`; i >= `0`; --i) {
897	CoinMemcpyN(rhs.index_ + start_[i], length_[i], index_ + start_[i]);
898	CoinMemcpyN(rhs.element_ + start_[i], length_[i], element_ + start_[i]);
899	}
900	}
901	} else {
902	copyOf(rhs);
903	}
904	}
905
906	//#############################################################################
907
908	// This method is essentially the same as minorAppendOrthoOrdered(). However,
909	// since we start from an empty matrix, lots of fluff can be avoided.
910
911	void
912	CoinPackedMatrix::reverseOrderedCopyOf(const CoinPackedMatrix& rhs)
913	{
914	if (this == &rhs) {
915	reverseOrdering();
916	return;
917	}
918
919	int i;
920	colOrdered_ = !rhs.colOrdered_;
921	majorDim_ = rhs.minorDim_;
922	minorDim_ = rhs.majorDim_;
923	size_ = rhs.size_;
924
925	if (size_ == `0`) {
926	// we still need to allocate starts and lengths
927	maxMajorDim_=majorDim_;
928	delete[] start_;
929	delete[] length_;
930	delete[] index_;
931	delete[] element_;
932	start_ = new CoinBigIndex[maxMajorDim_ + `1`];
933	length_ = new int[maxMajorDim_];
934	for (i = `0`; i < majorDim_; ++i) {
935	start_[i] = `0`;
936	length_[i]=`0`;
937	}
938	start_[majorDim_]=`0`;
939	index_ = new int[maxSize_];
940	element_ = new double[maxSize_];
941	return;
942	}
943
944
945	// Allocate sufficient space (resizeForAddingMinorVectors())
946
947	const int newMaxMajorDim_ =
948	CoinMax(maxMajorDim_, CoinLengthWithExtra(majorDim_, extraMajor_));
949
950	if (newMaxMajorDim_ > maxMajorDim_) {
951	maxMajorDim_ = newMaxMajorDim_;
952	delete[] start_;
953	delete[] length_;
954	start_ = new CoinBigIndex[maxMajorDim_ + `1`];
955	length_ = new int[maxMajorDim_];
956	}
957	// first compute how long each major-dimension vector will be
958	int * COIN_RESTRICT orthoLength = length_;
959	rhs.countOrthoLength(orthoLength);
960
961	start_[`0`] = `0`;
962	if (extraGap_ == `0`) {
963	for (i = `0`; i < majorDim_; ++i)
964	start_[i+`1`] = start_[i] + orthoLength[i];
965	} else {
966	const double eg = extraGap_;
967	for (i = `0`; i < majorDim_; ++i)
968	start_[i+`1`] = start_[i] + CoinLengthWithExtra(orthoLength[i], eg);
969	}
970
971	const CoinBigIndex newMaxSize =
972	CoinMax(maxSize_, CoinLengthWithExtra(getLastStart(), extraMajor_));
973
974	if (newMaxSize > maxSize_) {
975	maxSize_ = newMaxSize;
976	delete[] index_;
977	delete[] element_;
978	index_ = new int[maxSize_];
979	element_ = new double[maxSize_];
980	# ifdef ZEROFAULT
981	memset(index_,`0`,(maxSize_*sizeof(int))) ;
982	memset(element_,`0`,(maxSize_*sizeof(double))) ;
983	# endif
984	}
985
986	// now insert the entries of matrix
987
988	minorDim_ = rhs.majorDim_;
989	const CoinBigIndex * COIN_RESTRICT start = rhs.start_;
990	const int * COIN_RESTRICT index = rhs.index_;
991	const int * COIN_RESTRICT length = rhs.length_;
992	const double * COIN_RESTRICT element = rhs.element_;
993	assert (start[`0`]==`0`);
994	CoinBigIndex first = `0`;
995	for (i = `0`; i < minorDim_; ++i) {
996	CoinBigIndex last = first + length[i];
997	CoinBigIndex j = first;
998	first = start[i+`1`];
999	#if 0
1000	if (((last-j)&`1`)!=`0`) {
1001	const int ind = index[j];
1002	CoinBigIndex put = start_[ind];
1003	start_[ind] = put +`1`;
1004	element_[put] = element[j];
1005	index_[put] = i;
1006	j++;
1007	}
1008	for (; j != last; j+=`2`) {
1009	const int ind0 = index[j];
1010	CoinBigIndex put0 = start_[ind0];
1011	double value0=element[j];
1012	const int ind1 = index[j+`1`];
1013	CoinBigIndex put1 = start_[ind1];
1014	double value1=element[j+`1`];
1015	start_[ind0] = put0 +`1`;
1016	start_[ind1] = put1 +`1`;
1017	element_[put0] = value0;
1018	index_[put0] = i;
1019	element_[put1] = value1;
1020	index_[put1] = i;
1021	}
1022	#else
1023	for (; j != last; ++j) {
1024	const int ind = index[j];
1025	CoinBigIndex put = start_[ind];
1026	start_[ind] = put +`1`;
1027	element_[put] = element[j];
1028	index_[put] = i;
1029	}
1030	#endif
1031	}
1032	// and re-adjust start_
1033	for (i = `0`; i < majorDim_; ++i) {
1034	start_[i] -= length_[i];
1035	}
1036	}
1037
1038	//#############################################################################
1039
1040	void
1041	CoinPackedMatrix::assignMatrix(const bool colordered,
1042	const int minor, const int major,
1043	const CoinBigIndex numels,
1044	double & elem, int* *& ind,
1045	CoinBigIndex & start, int* *& len,
1046	const int maxmajor, const CoinBigIndex maxsize)
1047	{
1048	gutsOfDestructor();
1049	colOrdered_ = colordered;
1050	element_ = elem;
1051	index_ = ind;
1052	start_ = start;
1053	majorDim_ = major;
1054	minorDim_ = minor;
1055	size_ = numels;
1056	maxMajorDim_ = maxmajor != -`1` ? maxmajor : major;
1057	maxSize_ = maxsize != -`1` ? maxsize : numels;
1058	if (len == NULL) {
1059	delete [] length_;
1060	length_ = new int[maxMajorDim_];
1061	std::adjacent_difference(start + `1`, start + (major + `1`), length_);
1062	length_[`0`] -= start[`0`];
1063	} else {
1064	length_ = len;
1065	}
1066	elem = NULL;
1067	ind = NULL;
1068	start = NULL;
1069	len = NULL;
1070	}
1071
1072	//#############################################################################
1073
1074	CoinPackedMatrix &
1075	CoinPackedMatrix::operator=(const CoinPackedMatrix& rhs)
1076	{
1077	if (this != &rhs) {
1078	gutsOfDestructor();
1079	extraGap_=rhs.extraGap_;
1080	extraMajor_=rhs.extraMajor_;
1081	gutsOfOpEqual(rhs.colOrdered_,
1082	rhs.minorDim_, rhs.majorDim_, rhs.size_,
1083	rhs.element_, rhs.index_, rhs.start_, rhs.length_);
1084	}
1085	return *this;
1086	}
1087
1088	//#############################################################################
1089
1090	void
1091	CoinPackedMatrix::reverseOrdering()
1092	{
1093	CoinPackedMatrix m;
1094	m.extraGap_ = extraMajor_;
1095	m.extraMajor_ = extraGap_;
1096	m.reverseOrderedCopyOf(*this);
1097	swap(m);
1098	}
1099
1100	//-----------------------------------------------------------------------------
1101
1102	void
1103	CoinPackedMatrix::transpose()
1104	{
1105	colOrdered_ = ! colOrdered_;
1106	}
1107
1108	//-----------------------------------------------------------------------------
1109
1110	void
1111	CoinPackedMatrix::swap(CoinPackedMatrix& m)
1112	{
1113	std::swap(colOrdered_, m.colOrdered_);
1114	std::swap(extraGap_, m.extraGap_);
1115	std::swap(extraMajor_, m.extraMajor_);
1116	std::swap(element_, m.element_);
1117	std::swap(index_, m.index_);
1118	std::swap(start_, m.start_);
1119	std::swap(length_, m.length_);
1120	std::swap(majorDim_, m.majorDim_);
1121	std::swap(minorDim_, m.minorDim_);
1122	std::swap(size_, m.size_);
1123	std::swap(maxMajorDim_, m.maxMajorDim_);
1124	std::swap(maxSize_, m.maxSize_);
1125	}
1126
1127	//#############################################################################
1128	//#############################################################################
1129
1130	void
1131	CoinPackedMatrix::times(const double * x, double * y) const
1132	{
1133	if (colOrdered_)
1134	timesMajor(x, y);
1135	else
1136	timesMinor(x, y);
1137	}
1138
1139	//-----------------------------------------------------------------------------
1140	#ifndef CLP_NO_VECTOR
1141	void
1142	CoinPackedMatrix::times(const CoinPackedVectorBase& x, double * y) const
1143	{
1144	if (colOrdered_)
1145	timesMajor(x, y);
1146	else
1147	timesMinor(x, y);
1148	}
1149	#endif
1150	//-----------------------------------------------------------------------------
1151
1152	void
1153	CoinPackedMatrix::transposeTimes(const double * x, double * y) const
1154	{
1155	if (colOrdered_)
1156	timesMinor(x, y);
1157	else
1158	timesMajor(x, y);
1159	}
1160
1161	//-----------------------------------------------------------------------------
1162	#ifndef CLP_NO_VECTOR
1163	void
1164	CoinPackedMatrix::transposeTimes(const CoinPackedVectorBase& x, double * y) const
1165	{
1166	if (colOrdered_)
1167	timesMinor(x, y);
1168	else
1169	timesMajor(x, y);
1170	}
1171	#endif
1172	//#############################################################################
1173	//#############################################################################
1174	/ Count the number of entries in every minor-dimension vector and*
1175	fill in an array containing these lengths. /*
1176	void
1177	CoinPackedMatrix::countOrthoLength(int * orthoLength) const
1178	{
1179	CoinZeroN(orthoLength, minorDim_);
1180	if (size_!=start_[majorDim_]) {
1181	// has gaps
1182	for (int i = `0`; i <majorDim_ ; ++i) {
1183	const CoinBigIndex first = start_[i];
1184	const CoinBigIndex last = first + length_[i];
1185	for (CoinBigIndex j = first; j < last; ++j) {
1186	assert( index_[j] < minorDim_ && index_[j]>=`0`);
1187	++orthoLength[index_[j]];
1188	}
1189	}
1190	} else {
1191	// no gaps
1192	const CoinBigIndex last = start_[majorDim_];
1193	for (CoinBigIndex j = `0`; j < last; ++j) {
1194	assert( index_[j] < minorDim_ && index_[j]>=`0`);
1195	++orthoLength[index_[j]];
1196	}
1197	}
1198	}
1199
1200	int *
1201	CoinPackedMatrix::countOrthoLength() const
1202	{
1203	int * orthoLength = new int[minorDim_];
1204	countOrthoLength(orthoLength);
1205	return orthoLength;
1206	}
1207
1208	//#############################################################################
1209	/ Returns an array containing major indices. The array is*
1210	getNumElements long and if getVectorStarts() is 0,2,5 then
1211	the array would start 0,0,1,1,1,2...
1212	This method is provided to go back from a packed format
1213	to a triple format.
1214	The returned array is allocated with <code>new int[]</code>,
1215	free it with <code>delete[]</code>. /*
1216	int *
1217	CoinPackedMatrix::getMajorIndices() const
1218	{
1219	// Check valid
1220	if (!majorDim_\|\|start_[majorDim_]!=size_)
1221	return NULL;
1222	int * array = new int [size_];
1223	for (int i=`0`;i<majorDim_;i++) {
1224	for (CoinBigIndex k=start_[i];k<start_[i+`1`];k++)
1225	array[k]=i;
1226	}
1227	return array;
1228	}
1229	//#############################################################################
1230
1231	void
1232	CoinPackedMatrix::appendMajorVector(const int vecsize,
1233	const int *vecind,
1234	const double *vecelem)
1235	{
1236	#ifdef COIN_DEBUG
1237	for (int i = `0`; i < vecsize; ++i) {
1238	if (vecind[i] < `0` )
1239	throw CoinError("out of range index",
1240	"appendMajorVector", "CoinPackedMatrix");
1241	}
1242	#if 0
1243	if (std::find_if(vecind, vecind + vecsize,
1244	compose2(logical_or<bool>(),
1245	bind2nd(less<int>(), `0`),
1246	bind2nd(greater_equal<int>(), minorDim_))) !=
1247	vecind + vecsize)
1248	throw CoinError("out of range index",
1249	"appendMajorVector", "CoinPackedMatrix");
1250	#endif
1251	#endif
1252
1253	if (majorDim_ == maxMajorDim_ \|\| vecsize > maxSize_ - getLastStart()) {
1254	resizeForAddingMajorVectors(`1`, &vecsize);
1255	}
1256
1257	// got to get this again since it might change!
1258	const CoinBigIndex last = getLastStart();
1259
1260	// OK, now just append the major-dimension vector to the end
1261
1262	length_[majorDim_] = vecsize;
1263	CoinMemcpyN(vecind, vecsize, index_ + last);
1264	CoinMemcpyN(vecelem, vecsize, element_ + last);
1265	if (majorDim_ == `0`)
1266	start_[`0`] = `0`;
1267	start_[majorDim_ + `1`] =
1268	CoinMin(last + CoinLengthWithExtra(vecsize, extraGap_), maxSize_ );
1269
1270	// LL: Do we want to allow appending a vector that has more entries than
1271	// the current size?
1272	if (vecsize > `0`) {
1273	minorDim_ = CoinMax(minorDim_,
1274	(*std::max_element(vecind, vecind+vecsize)) + `1`);
1275	}
1276
1277	++majorDim_;
1278	size_ += vecsize;
1279	}
1280
1281	//-----------------------------------------------------------------------------
1282	#ifndef CLP_NO_VECTOR
1283	void
1284	CoinPackedMatrix::appendMajorVector(const CoinPackedVectorBase& vec)
1285	{
1286	appendMajorVector(vec.getNumElements(),
1287	vec.getIndices(), vec.getElements());
1288	}
1289	//-----------------------------------------------------------------------------
1290
1291	void
1292	CoinPackedMatrix::appendMajorVectors(const int numvecs,
1293	const CoinPackedVectorBase * const * vecs)
1294	{
1295	int i;
1296	CoinBigIndex nz = `0`;
1297	for (i = `0`; i < numvecs; ++i)
1298	nz += CoinLengthWithExtra(vecs[i]->getNumElements(), extraGap_);
1299	reserve(majorDim_ + numvecs, getLastStart() + nz);
1300	for (i = `0`; i < numvecs; ++i)
1301	appendMajorVector(*vecs[i]);
1302	}
1303	#endif
1304
1305	//#############################################################################
1306
1307	void
1308	CoinPackedMatrix::appendMinorVector(const int vecsize,
1309	const int *vecind,
1310	const double *vecelem)
1311	{
1312	if (vecsize == `0`) {
1313	++minorDim_; // empty row/column - still need to increase
1314	return;
1315	}
1316
1317	int i;
1318	#if COIN_COINUTILS_CHECKLEVEL > 3
1319	// Test if any of the indices are out of range
1320	for (i = `0`; i < vecsize; ++i) {
1321	if (vecind[i] < `0` \|\| vecind[i] >= majorDim_)
1322	throw CoinError("out of range index",
1323	"appendMinorVector", "CoinPackedMatrix");
1324	}
1325	// Test if there are duplicate indices
1326	int* sortedind = CoinCopyOfArray(vecind, vecsize);
1327	std::sort(sortedind, sortedind+vecsize);
1328	if (std::adjacent_find(sortedind, sortedind+vecsize) != sortedind+vecsize) {
1329	throw CoinError("identical indices",
1330	"appendMinorVector", "CoinPackedMatrix");
1331	}
1332	#endif
1333
1334	// test that there's a gap at the end of every major-dimension vector where
1335	// we want to add a new entry
1336
1337	for (i = vecsize - `1`; i >= `0`; --i) {
1338	const int j = vecind[i];
1339	if (start_[j] + length_[j] == start_[j+`1`])
1340	break;
1341	}
1342
1343	if (i >= `0`) {
1344	int * addedEntries = new int[majorDim_];
1345	memset(addedEntries, `0`, majorDim_ * sizeof(int));
1346	for (i = vecsize - `1`; i >= `0`; --i)
1347	addedEntries[vecind[i]] = `1`;
1348	resizeForAddingMinorVectors(addedEntries);
1349	delete[] addedEntries;
1350	}
1351
1352	// OK, now insert the entries of the minor-dimension vector
1353	for (i = vecsize - `1`; i >= `0`; --i) {
1354	const int j = vecind[i];
1355	const CoinBigIndex posj = start_[j] + (length_[j]++);
1356	index_[posj] = minorDim_;
1357	element_[posj] = vecelem[i];
1358	}
1359
1360	++minorDim_;
1361	size_ += vecsize;
1362	}
1363
1364	//-----------------------------------------------------------------------------
1365	#ifndef CLP_NO_VECTOR
1366	void
1367	CoinPackedMatrix::appendMinorVector(const CoinPackedVectorBase& vec)
1368	{
1369	appendMinorVector(vec.getNumElements(),
1370	vec.getIndices(), vec.getElements());
1371	}
1372
1373	//-----------------------------------------------------------------------------
1374
1375	void
1376	CoinPackedMatrix::appendMinorVectors(const int numvecs,
1377	const CoinPackedVectorBase * const * vecs)
1378	{
1379	if (numvecs == `0`)
1380	return;
1381
1382	int i;
1383
1384	int * addedEntries = new int[majorDim_];
1385	CoinZeroN(addedEntries, majorDim_);
1386	for (i = numvecs - `1`; i >= `0`; --i) {
1387	const int vecsize = vecs[i]->getNumElements();
1388	const int* vecind = vecs[i]->getIndices();
1389	for (int j = vecsize - `1`; j >= `0`; --j) {
1390	#ifdef COIN_DEBUG
1391	if (vecind[j] < `0` \|\| vecind[j] >= majorDim_)
1392	throw CoinError("out of range index", "appendMinorVectors",
1393	"CoinPackedMatrix");
1394	#endif
1395	++addedEntries[vecind[j]];
1396	}
1397	}
1398
1399	for (i = majorDim_ - `1`; i >= `0`; --i) {
1400	if (start_[i] + length_[i] + addedEntries[i] > start_[i+`1`])
1401	break;
1402	}
1403	if (i >= `0`)
1404	resizeForAddingMinorVectors(addedEntries);
1405	delete[] addedEntries;
1406
1407	// now insert the entries of the vectors
1408	for (i = `0`; i < numvecs; ++i) {
1409	const int vecsize = vecs[i]->getNumElements();
1410	const int* vecind = vecs[i]->getIndices();
1411	const double* vecelem = vecs[i]->getElements();
1412	for (int j = vecsize - `1`; j >= `0`; --j) {
1413	const int ind = vecind[j];
1414	element_[start_[ind] + length_[ind]] = vecelem[j];
1415	index_[start_[ind] + (length_[ind]++)] = minorDim_;
1416	}
1417	++minorDim_;
1418	size_ += vecsize;
1419	}
1420	}
1421	#endif
1422
1423	//#############################################################################
1424	//#############################################################################
1425
1426	void
1427	CoinPackedMatrix::majorAppendSameOrdered(const CoinPackedMatrix& matrix)
1428	{
1429	if (minorDim_ != matrix.minorDim_) {
1430	throw CoinError ("dimension mismatch", "rightAppendSameOrdered",
1431	"CoinPackedMatrix");
1432	}
1433	if (matrix.majorDim_ == `0`)
1434	return;
1435
1436	int i;
1437	if (majorDim_ + matrix.majorDim_ > maxMajorDim_ \|\|
1438	getLastStart() + matrix.getLastStart() > maxSize_) {
1439	// we got to resize before we add. note that the resizing method
1440	// properly fills out start_ and length_ for the major-dimension
1441	// vectors to be added!
1442	resizeForAddingMajorVectors(matrix.majorDim_, matrix.length_);
1443	start_ += majorDim_;
1444	for (i = `0`; i < matrix.majorDim_; ++i) {
1445	const int l = matrix.length_[i];
1446	CoinMemcpyN(matrix.index_ + matrix.start_[i], l,
1447	index_ + start_[i]);
1448	CoinMemcpyN(matrix.element_ + matrix.start_[i], l,
1449	element_ + start_[i]);
1450	}
1451	start_ -= majorDim_;
1452	} else {
1453	start_ += majorDim_;
1454	length_ += majorDim_;
1455	for (i = `0`; i < matrix.majorDim_; ++i) {
1456	const int l = matrix.length_[i];
1457	CoinMemcpyN(matrix.index_ + matrix.start_[i], l,
1458	index_ + start_[i]);
1459	CoinMemcpyN(matrix.element_ + matrix.start_[i], l,
1460	element_ + start_[i]);
1461	start_[i+`1`] = start_[i] + matrix.start_[i+`1`] - matrix.start_[i];
1462	length_[i] = l;
1463	}
1464	start_ -= majorDim_;
1465	length_ -= majorDim_;
1466	}
1467	majorDim_ += matrix.majorDim_;
1468	size_ += matrix.size_;
1469	}
1470
1471	//-----------------------------------------------------------------------------
1472
1473	void
1474	CoinPackedMatrix::minorAppendSameOrdered(const CoinPackedMatrix& matrix)
1475	{
1476	if (majorDim_ != matrix.majorDim_) {
1477	throw CoinError ("dimension mismatch", "bottomAppendSameOrdered",
1478	"CoinPackedMatrix");
1479	}
1480	if (matrix.minorDim_ == `0`)
1481	return;
1482
1483	int i;
1484	for (i = majorDim_ - `1`; i >= `0`; --i) {
1485	if (start_[i] + length_[i] + matrix.length_[i] > start_[i+`1`])
1486	break;
1487	}
1488	if (i >= `0`)
1489	resizeForAddingMinorVectors(matrix.length_);
1490
1491	// now insert the entries of matrix
1492	for (i = majorDim_ - `1`; i >= `0`; --i) {
1493	const int l = matrix.length_[i];
1494	std::transform(matrix.index_ + matrix.start_[i],
1495	matrix.index_ + (matrix.start_[i] + l),
1496	index_ + (start_[i] + length_[i]),
1497	std::bind2nd(std::plus<int>(), minorDim_));
1498	CoinMemcpyN(matrix.element_ + matrix.start_[i], l,
1499	element_ + (start_[i] + length_[i]));
1500	length_[i] += l;
1501	}
1502	minorDim_ += matrix.minorDim_;
1503	size_ += matrix.size_;
1504	}
1505
1506	//-----------------------------------------------------------------------------
1507
1508	void
1509	CoinPackedMatrix::majorAppendOrthoOrdered(const CoinPackedMatrix& matrix)
1510	{
1511	if (minorDim_ != matrix.majorDim_) {
1512	throw CoinError ("dimension mismatch", "majorAppendOrthoOrdered",
1513	"CoinPackedMatrix");
1514	}
1515	if (matrix.majorDim_ == `0`)
1516	return;
1517
1518	int i;
1519	CoinBigIndex j;
1520	// this trickery is needed because MSVC++ is not willing to delete[] a
1521	// 'const int '*
1522	int * orthoLengthPtr = matrix.countOrthoLength();
1523	const int * orthoLength = orthoLengthPtr;
1524
1525	if (majorDim_ + matrix.minorDim_ > maxMajorDim_) {
1526	resizeForAddingMajorVectors(matrix.minorDim_, orthoLength);
1527	} else {
1528	const double extra_gap = extraGap_;
1529	start_ += majorDim_;
1530	for (i = `0`; i < matrix.minorDim_ ; ++i) {
1531	start_[i+`1`] = start_[i] + CoinLengthWithExtra(orthoLength[i], extra_gap);
1532	}
1533	start_ -= majorDim_;
1534	if (start_[majorDim_ + matrix.minorDim_] > maxSize_) {
1535	resizeForAddingMajorVectors(matrix.minorDim_, orthoLength);
1536	}
1537	}
1538	// At this point everything is big enough to accommodate the new entries.
1539	// Also, start_ is set to the correct starting points for all the new
1540	// major-dimension vectors. The length of the new major-dimension vectors
1541	// may or may not be correctly set. Hence we just zero them out and they'll
1542	// be set when the entries are actually added below.
1543
1544	start_ += majorDim_;
1545	length_ += majorDim_;
1546
1547	CoinZeroN(length_, matrix.minorDim_);
1548
1549	for (i = `0`; i < matrix.majorDim_; ++i) {
1550	const CoinBigIndex last = matrix.getVectorLast(i);
1551	for (j = matrix.getVectorFirst(i); j < last; ++j) {
1552	const int ind = matrix.index_[j];
1553	element_[start_[ind] + length_[ind]] = matrix.element_[j];
1554	index_[start_[ind] + (length_[ind]++)] = i;
1555	}
1556	}
1557
1558	length_ -= majorDim_;
1559	start_ -= majorDim_;
1560
1561	// We need to update majorDim_ and size_. We can just add in from matrix
1562	majorDim_ += matrix.minorDim_;
1563	size_ += matrix.size_;
1564
1565	delete[] orthoLengthPtr;
1566	}
1567
1568	//-----------------------------------------------------------------------------
1569
1570	void
1571	CoinPackedMatrix::minorAppendOrthoOrdered(const CoinPackedMatrix& matrix)
1572	{
1573	if (majorDim_ != matrix.minorDim_) {
1574	throw CoinError ("dimension mismatch", "bottomAppendOrthoOrdered",
1575	"CoinPackedMatrix");
1576	}
1577	if (matrix.majorDim_ == `0`)
1578	return;
1579
1580	int i;
1581	// first compute how many entries will be added to each major-dimension
1582	// vector, and if needed, resize the matrix to accommodate all
1583	// this trickery is needed because MSVC++ is not willing to delete[] a
1584	// 'const int '*
1585	int * addedEntriesPtr = matrix.countOrthoLength();
1586	const int * addedEntries = addedEntriesPtr;
1587	for (i = majorDim_ - `1`; i >= `0`; --i) {
1588	if (start_[i] + length_[i] + addedEntries[i] > start_[i+`1`])
1589	break;
1590	}
1591	if (i >= `0`)
1592	resizeForAddingMinorVectors(addedEntries);
1593	delete[] addedEntriesPtr;
1594
1595	// now insert the entries of matrix
1596	for (i = `0`; i < matrix.majorDim_; ++i) {
1597	const CoinBigIndex last = matrix.getVectorLast(i);
1598	for (CoinBigIndex j = matrix.getVectorFirst(i); j != last; ++j) {
1599	const int ind = matrix.index_[j];
1600	element_[start_[ind] + length_[ind]] = matrix.element_[j];
1601	index_[start_[ind] + (length_[ind]++)] = minorDim_;
1602	}
1603	++minorDim_;
1604	}
1605	size_ += matrix.size_;
1606	}
1607
1608	//#############################################################################
1609	//#############################################################################
1610
1611	void
1612	CoinPackedMatrix::deleteMajorVectors(const int numDel,
1613	const int * indDel)
1614	{
1615	if (numDel == majorDim_) {
1616	// everything is deleted
1617	majorDim_ = `0`;
1618	minorDim_ = `0`;
1619	size_ = `0`;
1620	// Get rid of memory as well
1621	maxMajorDim_ = `0`;
1622	delete [] length_;
1623	length_ = NULL;
1624	delete [] start_;
1625	start_ = new CoinBigIndex[`1`];
1626	start_[`0`]=`0`;
1627	delete [] element_;
1628	element_=NULL;
1629	delete [] index_;
1630	index_=NULL;
1631	maxSize_ = `0`;
1632	return;
1633	}
1634
1635	if (!extraGap_&&!extraMajor_) {
1636	// See if this is faster
1637	char * keep = new char[majorDim_];
1638	memset(keep,`1`,majorDim_);
1639	for (int i=`0`;i<numDel;i++) {
1640	int k=indDel[i];
1641	assert (k>=`0`&&k<majorDim_&&keep[k]);
1642	keep[k]=`0`;
1643	}
1644	int n;
1645	// find first
1646	for (n=`0`;n<majorDim_;n++) {
1647	if (!keep[n])
1648	break;
1649	}
1650	size_=start_[n];
1651	for (int i=n;i<majorDim_;i++) {
1652	if (keep[i]) {
1653	int length = length_[i];
1654	length_[n]=length;
1655	for (CoinBigIndex j=start_[i];j<start_[i+`1`];j++) {
1656	element_[size_]=element_[j];
1657	index_[size_++]=index_[j];
1658	}
1659	start_[++n]=size_;
1660	}
1661	}
1662	majorDim_=n;
1663	delete [] keep;
1664	} else {
1665	int *sortedDelPtr = CoinTestIndexSet(numDel, indDel, majorDim_,
1666	"deleteMajorVectors");
1667	const int * sortedDel = sortedDelPtr == `0` ? indDel : sortedDelPtr;
1668
1669	CoinBigIndex deleted = `0`;
1670	const int last = numDel - `1`;
1671	for (int i = `0`; i < last; ++i) {
1672	const int ind = sortedDel[i];
1673	const int ind1 = sortedDel[i+`1`];
1674	deleted += length_[ind];
1675	if (ind1 - ind > `1`) {
1676	CoinCopy(start_ + (ind + `1`), start_ + ind1, start_ + (ind - i));
1677	CoinCopy(length_ + (ind + `1`), length_ + ind1, length_ + (ind - i));
1678	}
1679	}
1680
1681	// copy the last block of length_ and start_
1682	const int ind = sortedDel[last];
1683	deleted += length_[ind];
1684	if (sortedDel[last] != majorDim_ - `1`) {
1685	const int ind1 = majorDim_;
1686	CoinCopy(start_ + (ind + `1`), start_ + ind1, start_ + (ind - last));
1687	CoinCopy(length_ + (ind + `1`), length_ + ind1, length_ + (ind - last));
1688	}
1689	majorDim_ -= numDel;
1690	const int lastlength = CoinLengthWithExtra(length_[majorDim_-`1`], extraGap_);
1691	start_[majorDim_] = CoinMin(start_[majorDim_-`1`] + lastlength, maxSize_);
1692	size_ -= deleted;
1693
1694	// if the very first major vector was deleted then copy the new first major
1695	// vector to the beginning to make certain that start_[0] is 0. This may
1696	// not be necessary, but better safe than sorry...
1697	if (sortedDel[`0`] == `0`) {
1698	CoinCopyN(index_ + start_[`0`], length_[`0`], index_);
1699	CoinCopyN(element_ + start_[`0`], length_[`0`], element_);
1700	start_[`0`] = `0`;
1701	}
1702
1703	delete[] sortedDelPtr;
1704	}
1705	}
1706
1707	//#############################################################################
1708
1709	void
1710	CoinPackedMatrix::deleteMinorVectors(const int numDel,
1711	const int * indDel)
1712	{
1713	if (numDel == minorDim_) {
1714	// everything is deleted
1715	minorDim_ = `0`;
1716	size_ = `0`;
1717	// Get rid of as much memory as possible
1718	memset(length_,`0`,majorDim_*sizeof(int));
1719	memset(start_,`0`,(majorDim_+`1`)*sizeof(CoinBigIndex ));
1720	delete [] element_;
1721	element_=NULL;
1722	delete [] index_;
1723	index_=NULL;
1724	maxSize_ = `0`;
1725	return;
1726	}
1727	int i, j, k;
1728
1729	// first compute the new index of every row
1730	int* newindexPtr = new int[minorDim_];
1731	CoinZeroN(newindexPtr, minorDim_);
1732	for (j = `0`; j < numDel; ++j) {
1733	const int ind = indDel[j];
1734	#ifdef COIN_DEBUG
1735	if (ind < `0` \|\| ind >= minorDim_)
1736	throw CoinError("out of range index",
1737	"deleteMinorVectors", "CoinPackedMatrix");
1738	if (newindexPtr[ind] == -`1`)
1739	throw CoinError("duplicate index",
1740	"deleteMinorVectors", "CoinPackedMatrix");
1741	#endif
1742	newindexPtr[ind] = -`1`;
1743	}
1744	for (i = `0`, k = `0`; i < minorDim_; ++i) {
1745	if (newindexPtr[i] != -`1`) {
1746	newindexPtr[i] = k++;
1747	}
1748	}
1749	// Now crawl through the matrix
1750	const int * newindex = newindexPtr;
1751	#ifdef TAKEOUT
1752	int mcount[`400`];
1753	memset(mcount,`0`,`400`*sizeof(int));
1754	for (i = `0`; i < majorDim_; ++i) {
1755	int * index = index_ + start_[i];
1756	double * elem = element_ + start_[i];
1757	const int length_i = length_[i];
1758	for (j = `0`, k = `0`; j < length_i; ++j) {
1759	mcount[index[j]]++;
1760	}
1761	}
1762	for (i=`0`;i<minorDim_;i++) {
1763	if (mcount[i]==`10`\|\|mcount[i]==`15`) {
1764	if (newindex[i]>=`0`)
1765	printf("Keeping original row %d (new %d) with count of %d\n",
1766	i,newindex[i],mcount[i]);
1767	else
1768	printf("deleting row %d with count of %d\n",
1769	i,mcount[i]);
1770	}
1771	}
1772	#endif
1773	if (!extraGap_) {
1774	// pack down
1775	size_=`0`;
1776	for (i = `0`; i < majorDim_; ++i) {
1777	int * index = index_ + start_[i];
1778	double * elem = element_ + start_[i];
1779	start_[i]=size_;
1780	const int length_i = length_[i];
1781	for (j = `0`; j < length_i; ++j) {
1782	const int ind = newindex[index[j]];
1783	if (ind >= `0`) {
1784	index_[size_] = ind;
1785	element_[size_++] = elem[j];
1786	}
1787	}
1788	length_[i] = size_-start_[i];
1789	}
1790	start_[majorDim_]=size_;
1791	} else {
1792	int deleted = `0`;
1793	for (i = `0`; i < majorDim_; ++i) {
1794	int * index = index_ + start_[i];
1795	double * elem = element_ + start_[i];
1796	const int length_i = length_[i];
1797	for (j = `0`, k = `0`; j < length_i; ++j) {
1798	const int ind = newindex[index[j]];
1799	if (ind != -`1`) {
1800	index[k] = ind;
1801	elem[k++] = elem[j];
1802	}
1803	}
1804	deleted += length_i - k;
1805	length_[i] = k;
1806	}
1807	size_ -= deleted;
1808	}
1809
1810	delete[] newindexPtr;
1811
1812	minorDim_ -= numDel;
1813	}
1814
1815	//#############################################################################
1816	//#############################################################################
1817
1818	void
1819	CoinPackedMatrix::timesMajor(const double * x, double * y) const
1820	{
1821	memset(y, `0`, minorDim_ * sizeof(double));
1822	for (int i = majorDim_ - `1`; i >= `0`; --i) {
1823	const double x_i = x[i];
1824	if (x_i != `0.0`) {
1825	const CoinBigIndex last = getVectorLast(i);
1826	for (CoinBigIndex j = getVectorFirst(i); j < last; ++j)
1827	y[index_[j]] += x_i * element_[j];
1828	}
1829	}
1830	}
1831
1832	//-----------------------------------------------------------------------------
1833	#ifndef CLP_NO_VECTOR
1834	void
1835	CoinPackedMatrix::timesMajor(const CoinPackedVectorBase& x, double * y) const
1836	{
1837	memset(y, `0`, minorDim_ * sizeof(double));
1838	for (CoinBigIndex i = x.getNumElements() - `1`; i >= `0`; --i) {
1839	const double x_i = x.getElements()[i];
1840	if (x_i != `0.0`) {
1841	const int ind = x.getIndices()[i];
1842	const CoinBigIndex last = getVectorLast(ind);
1843	for (CoinBigIndex j = getVectorFirst(ind); j < last; ++j)
1844	y[index_[j]] += x_i * element_[j];
1845	}
1846	}
1847	}
1848	#endif
1849	//-----------------------------------------------------------------------------
1850
1851	void
1852	CoinPackedMatrix::timesMinor(const double * x, double * y) const
1853	{
1854	memset(y, `0`, majorDim_ * sizeof(double));
1855	for (int i = majorDim_ - `1`; i >= `0`; --i) {
1856	double y_i = `0`;
1857	const CoinBigIndex last = getVectorLast(i);
1858	for (CoinBigIndex j = getVectorFirst(i); j < last; ++j)
1859	y_i += x[index_[j]] * element_[j];
1860	y[i] = y_i;
1861	}
1862	}
1863
1864	//-----------------------------------------------------------------------------
1865	#ifndef CLP_NO_VECTOR
1866	void
1867	CoinPackedMatrix::timesMinor(const CoinPackedVectorBase& x, double * y) const
1868	{
1869	memset(y, `0`, majorDim_ * sizeof(double));
1870	for (int i = majorDim_ - `1`; i >= `0`; --i) {
1871	double y_i = `0`;
1872	const CoinBigIndex last = getVectorLast(i);
1873	for (CoinBigIndex j = getVectorFirst(i); j < last; ++j)
1874	y_i += x [index_[j]] * element_[j];
1875	y[i] = y_i;
1876	}
1877	}
1878	#endif
1879	//#############################################################################
1880	//#############################################################################
1881
1882	CoinPackedMatrix::CoinPackedMatrix() :
1883	colOrdered_(true),
1884	extraGap_(`0.0`),
1885	extraMajor_(`0.0`),
1886	element_(`0`),
1887	index_(`0`),
1888	length_(`0`),
1889	majorDim_(`0`),
1890	minorDim_(`0`),
1891	size_(`0`),
1892	maxMajorDim_(`0`),
1893	maxSize_(`0`)
1894	{
1895	start_ = new CoinBigIndex[`1`];
1896	start_[`0`] = `0`;
1897	}
1898
1899	//-----------------------------------------------------------------------------
1900
1901	CoinPackedMatrix::CoinPackedMatrix(const bool colordered,
1902	const double extraMajor,
1903	const double extraGap) :
1904	colOrdered_(colordered),
1905	extraGap_(extraGap),
1906	extraMajor_(extraMajor),
1907	element_(`0`),
1908	index_(`0`),
1909	length_(`0`),
1910	majorDim_(`0`),
1911	minorDim_(`0`),
1912	size_(`0`),
1913	maxMajorDim_(`0`),
1914	maxSize_(`0`)
1915	{
1916	start_ = new CoinBigIndex[`1`];
1917	start_[`0`] = `0`;
1918	}
1919
1920	//-----------------------------------------------------------------------------
1921
1922	CoinPackedMatrix::CoinPackedMatrix(const bool colordered,
1923	const int minor, const int major,
1924	const CoinBigIndex numels,
1925	const double * elem, const int * ind,
1926	const CoinBigIndex * start, const int * len,
1927	const double extraMajor,
1928	const double extraGap) :
1929	colOrdered_(colordered),
1930	extraGap_(extraGap),
1931	extraMajor_(extraMajor),
1932	element_(NULL),
1933	index_(NULL),
1934	start_(NULL),
1935	length_(NULL),
1936	majorDim_(`0`),
1937	minorDim_(`0`),
1938	size_(`0`),
1939	maxMajorDim_(`0`),
1940	maxSize_(`0`)
1941	{
1942	gutsOfOpEqual(colordered, minor, major, numels, elem, ind, start, len);
1943	}
1944
1945	//-----------------------------------------------------------------------------
1946
1947	CoinPackedMatrix::CoinPackedMatrix(const bool colordered,
1948	const int minor, const int major,
1949	const CoinBigIndex numels,
1950	const double * elem, const int * ind,
1951	const CoinBigIndex * start, const int * len) :
1952	colOrdered_(colordered),
1953	extraGap_(`0.0`),
1954	extraMajor_(`0.0`),
1955	element_(NULL),
1956	index_(NULL),
1957	start_(NULL),
1958	length_(NULL),
1959	majorDim_(`0`),
1960	minorDim_(`0`),
1961	size_(`0`),
1962	maxMajorDim_(`0`),
1963	maxSize_(`0`)
1964	{
1965	gutsOfOpEqual(colordered, minor, major, numels, elem, ind, start, len);
1966	}
1967
1968	//-----------------------------------------------------------------------------
1969	// makes column ordered from triplets and takes out duplicates
1970	// will be sorted
1971	//
1972	// This is an interesting in-place sorting algorithm;
1973	// We have triples, and want to sort them so that triples with the same column
1974	// are adjacent.
1975	// We begin by computing how many entries there are for each column (columnCount)
1976	// and using that to compute where each set of column entries will end* (startColumn).*
1977	// As we drop entries into place, startColumn is decremented until it contains
1978	// the position where the column entries start.
1979	// The invalid column index -2 means there's a "hole" in that position;
1980	// the invalid column index -1 means the entry in that spot is "where it wants to go".
1981	// Initially, no one is where they want to go.
1982	// Going back to front,
1983	// if that entry is where it wants to go
1984	// then leave it there
1985	// otherwise pick it up (which leaves a hole), and
1986	// for as long as you have an entry in your right hand,
1987	// - pick up the entry (with your left hand) in the position where the one in
1988	// your right hand wants to go;
1989	// - pass the entry in your left hand to your right hand;
1990	// - was that entry really just the "hole"? If so, stop.
1991	// It could be that all the entries get shuffled in the first loop iteration
1992	// and all the rest just confirm that everyone is happy where they are.
1993	// We never move an entry that is where it wants to go, so entries are moved at
1994	// most once. They may not change position if they happen to initially be
1995	// where they want to go when the for loop gets to them.
1996	// It depends on how many subpermutations the triples initially defined.
1997	// Each while loop takes care of one permutation.
1998	// The while loop has to stop, because each time around we mark one entry as happy.
1999	// We can't run into a happy entry, because we are decrementing the startColumn
2000	// all the time, so we must be running into new entries.
2001	// Once we've processed all the slots for a column, it cannot be the case that
2002	// there are any others that want to go there.
2003	// This all means that we eventually must run into the hole.
2004	CoinPackedMatrix::CoinPackedMatrix(
2005	const bool colordered,
2006	const int * indexRow ,
2007	const int * indexColumn,
2008	const double * element,
2009	CoinBigIndex numberElements )
2010	:
2011	colOrdered_(colordered),
2012	extraGap_(`0.0`),
2013	extraMajor_(`0.0`),
2014	element_(NULL),
2015	index_(NULL),
2016	start_(NULL),
2017	length_(NULL),
2018	majorDim_(`0`),
2019	minorDim_(`0`),
2020	size_(`0`),
2021	maxMajorDim_(`0`),
2022	maxSize_(`0`)
2023	{
2024	CoinAbsFltEq eq;
2025	int * colIndices = new int[numberElements];
2026	int * rowIndices = new int[numberElements];
2027	double * elements = new double[numberElements];
2028	CoinCopyN(element,numberElements,elements);
2029	if ( colordered ) {
2030	CoinCopyN(indexColumn,numberElements,colIndices);
2031	CoinCopyN(indexRow,numberElements,rowIndices);
2032	}
2033	else {
2034	CoinCopyN(indexColumn,numberElements,rowIndices);
2035	CoinCopyN(indexRow,numberElements,colIndices);
2036	}
2037
2038	int numberRows;
2039	int numberColumns;
2040	if (numberElements ) {
2041	numberRows = *std::max_element(rowIndices,rowIndices+numberElements)+`1`;
2042	numberColumns = *std::max_element(colIndices,colIndices+numberElements)+`1`;
2043	} else {
2044	numberRows = `0`;
2045	numberColumns = `0`;
2046	}
2047	int * rowCount = new int[numberRows];
2048	int * columnCount = new int[numberColumns];
2049	CoinBigIndex * startColumn = new CoinBigIndex[numberColumns+`1`];
2050	int * lengths = new int[numberColumns+`1`];
2051
2052	int iColumn,i;
2053	CoinBigIndex k;
2054	for (i=`0`;i<numberRows;i++) {
2055	rowCount[i]=`0`;
2056	}
2057	for (i=`0`;i<numberColumns;i++) {
2058	columnCount[i]=`0`;
2059	}
2060	for (i=`0`;i<numberElements;i++) {
2061	int iRow=rowIndices[i];
2062	int iColumn=colIndices[i];
2063	rowCount[iRow]++;
2064	columnCount[iColumn]++;
2065	}
2066	CoinBigIndex iCount=`0`;
2067	for (iColumn=`0`;iColumn<numberColumns;iColumn++) {
2068	/ position after end of Column /
2069	iCount+=columnCount[iColumn];
2070	startColumn[iColumn]=iCount;
2071	} / endfor /
2072	startColumn[iColumn]=iCount;
2073	for (k=numberElements-`1`;k>=`0`;k--) {
2074	iColumn=colIndices[k];
2075	if (iColumn>=`0`) {
2076	/ pick up the entry with your right hand /
2077	double value = elements[k];
2078	int iRow=rowIndices[k];
2079	colIndices[k]=-`2`; / the hole /
2080
2081	while (`1`) {
2082	/ pick this up with your left /
2083	CoinBigIndex iLook=startColumn[iColumn]-`1`;
2084	double valueSave=elements[iLook];
2085	int iColumnSave=colIndices[iLook];
2086	int iRowSave=rowIndices[iLook];
2087
2088	/ put the right-hand entry where it wanted to go /
2089	startColumn[iColumn]=iLook;
2090	elements[iLook]=value;
2091	rowIndices[iLook]=iRow;
2092	colIndices[iLook]=-`1`; / mark it as being where it wants to be /
2093
2094	/ there was something there /
2095	if (iColumnSave>=`0`) {
2096	iColumn=iColumnSave;
2097	value=valueSave;
2098	iRow=iRowSave;
2099	} else if (iColumnSave == -`2`) { / that was the hole /
2100	break;
2101	} else {
2102	assert(`1`==`0`); / should never happen /
2103	}
2104	/ endif /
2105	} / endwhile /
2106	} / endif /
2107	} / endfor /
2108
2109	/ now pack the elements and combine entries with the same row and column /
2110	/ also, drop entries with "small" coefficients /
2111	numberElements=`0`;
2112	for (iColumn=`0`;iColumn<numberColumns;iColumn++) {
2113	CoinBigIndex start=startColumn[iColumn];
2114	CoinBigIndex end =startColumn[iColumn+`1`];
2115	lengths[iColumn]=`0`;
2116	startColumn[iColumn]=numberElements;
2117	if (end>start) {
2118	int lastRow;
2119	double lastValue;
2120	// sorts on indices dragging elements with
2121	CoinSort_2(rowIndices+start,rowIndices+end,elements+start,CoinFirstLess_2<int, double>());
2122	lastRow=rowIndices[start];
2123	lastValue=elements[start];
2124	for (i=start+`1`;i<end;i++) {
2125	int iRow=rowIndices[i];
2126	double value=elements[i];
2127	if (iRow>lastRow) {
2128	//if(fabs(lastValue)>tolerance) {
2129	if(!eq (lastValue,`0.0`)) {
2130	rowIndices[numberElements]=lastRow;
2131	elements[numberElements]=lastValue;
2132	numberElements++;
2133	lengths[iColumn]++;
2134	}
2135	lastRow=iRow;
2136	lastValue=value;
2137	} else {
2138	lastValue+=value;
2139	} / endif /
2140	} / endfor /
2141	//if(fabs(lastValue)>tolerance) {
2142	if(!eq (lastValue,`0.0`)) {
2143	rowIndices[numberElements]=lastRow;
2144	elements[numberElements]=lastValue;
2145	numberElements++;
2146	lengths[iColumn]++;
2147	}
2148	}
2149	} / endfor /
2150	startColumn[numberColumns]=numberElements;
2151	#if 0
2152	gutsOfOpEqual(colordered,numberRows,numberColumns,numberElements,elements,rowIndices,startColumn,lengths);
2153
2154	delete [] rowCount;
2155	delete [] columnCount;
2156	delete [] startColumn;
2157	delete [] lengths;
2158
2159	delete [] colIndices;
2160	delete [] rowIndices;
2161	delete [] elements;
2162	#else
2163	assignMatrix(colordered,numberRows,numberColumns,numberElements,
2164	elements,rowIndices,startColumn,lengths);
2165	delete [] rowCount;
2166	delete [] columnCount;
2167	delete [] lengths;
2168	delete [] colIndices;
2169	#endif
2170
2171	}
2172
2173	//-----------------------------------------------------------------------------
2174
2175	CoinPackedMatrix::CoinPackedMatrix (const CoinPackedMatrix & rhs) :
2176	colOrdered_(true),
2177	extraGap_(`0.0`),
2178	extraMajor_(`0.0`),
2179	element_(`0`),
2180	index_(`0`),
2181	start_(`0`),
2182	length_(`0`),
2183	majorDim_(`0`),
2184	minorDim_(`0`),
2185	size_(`0`),
2186	maxMajorDim_(`0`),
2187	maxSize_(`0`)
2188	{
2189	bool hasGaps = rhs.size_<rhs.start_[rhs.majorDim_];
2190	if (!hasGaps&&!rhs.extraMajor_) {
2191	gutsOfCopyOfNoGaps(rhs.colOrdered_,
2192	rhs.minorDim_, rhs.majorDim_,
2193	rhs.element_, rhs.index_, rhs.start_);
2194	} else {
2195	gutsOfCopyOf(rhs.colOrdered_,
2196	rhs.minorDim_, rhs.majorDim_, rhs.size_,
2197	rhs.element_, rhs.index_, rhs.start_, rhs.length_,
2198	rhs.extraMajor_, rhs.extraGap_);
2199	}
2200	}
2201	/ Copy constructor - fine tuning - allowing extra space and/or reverse ordering.*
2202	extraForMajor is exact extra after any possible reverse ordering.
2203	extraMajor_ and extraGap_ set to zero.
2204	*/
2205	CoinPackedMatrix::CoinPackedMatrix(const CoinPackedMatrix& rhs, int extraForMajor,
2206	int extraElements, bool reverseOrdering)
2207	: colOrdered_(rhs.colOrdered_),
2208	extraGap_(`0`),
2209	extraMajor_(`0`),
2210	element_(`0`),
2211	index_(`0`),
2212	start_(`0`),
2213	length_(`0`),
2214	majorDim_(rhs.majorDim_),
2215	minorDim_(rhs.minorDim_),
2216	size_(rhs.size_),
2217	maxMajorDim_(`0`),
2218	maxSize_(`0`)
2219	{
2220	if (!reverseOrdering) {
2221	if (extraForMajor>=`0`) {
2222	maxMajorDim_ = majorDim_+ extraForMajor;
2223	maxSize_ = size_ + extraElements;
2224	assert (maxMajorDim_>`0`);
2225	assert (maxSize_>`0`);
2226	length_ = new int[maxMajorDim_];
2227	CoinMemcpyN(rhs.length_, majorDim_, length_);
2228	start_ = new CoinBigIndex[maxMajorDim_+`1`];
2229	element_ = new double[maxSize_];
2230	index_ = new int[maxSize_];
2231	bool hasGaps = rhs.size_<rhs.start_[rhs.majorDim_];
2232	if (hasGaps) {
2233	// we can't just simply memcpy these content over, because that can
2234	// upset memory debuggers like purify if there were gaps and those gaps
2235	// were uninitialized memory blocks
2236	CoinBigIndex size=`0`;
2237	for (int i = `0` ; i < majorDim_ ; i++) {
2238	start_[i]=size;
2239	CoinMemcpyN(rhs.index_ + rhs.start_[i], length_[i], index_ + size);
2240	CoinMemcpyN(rhs.element_ + rhs.start_[i], length_[i], element_ + size);
2241	size += length_[i];
2242	}
2243	start_[majorDim_]=size;
2244	assert (size_==size);
2245	} else {
2246	CoinMemcpyN(rhs.start_, majorDim_+`1`, start_);
2247	CoinMemcpyN(rhs.index_, size_, index_);
2248	CoinMemcpyN(rhs.element_, size_, element_ );
2249	}
2250	} else {
2251	// take out small and gaps
2252	maxMajorDim_ = majorDim_;
2253	maxSize_ = size_;
2254	if (maxMajorDim_>`0`) {
2255	length_ = new int[maxMajorDim_];
2256	start_ = new CoinBigIndex[maxMajorDim_+`1`];
2257	if (maxSize_>`0`) {
2258	element_ = new double[maxSize_];
2259	index_ = new int[maxSize_];
2260	}
2261	CoinBigIndex size=`0`;
2262	const double * oldElement = rhs.element_;
2263	const CoinBigIndex * oldStart = rhs.start_;
2264	const int * oldIndex = rhs.index_;
2265	const int * oldLength = rhs.length_;
2266	CoinBigIndex tooSmallCount=`0`;
2267	for (int i = `0` ; i < majorDim_ ; i++) {
2268	start_[i]=size;
2269	for (CoinBigIndex j=oldStart[i];
2270	j<oldStart[i]+oldLength[i];j++) {
2271	double value = oldElement[j];
2272	if (fabs(value)>`1.0e-21`) {
2273	element_[size]=value;
2274	index_[size++]=oldIndex[j];
2275	} else {
2276	tooSmallCount++;
2277	}
2278	}
2279	length_[i]=size-start_[i];
2280	}
2281	start_[majorDim_]=size;
2282	assert (size_==size+tooSmallCount);
2283	size_ = size;
2284	} else {
2285	start_ = new CoinBigIndex[`1`];
2286	start_[`0`]=`0`;
2287	}
2288	}
2289	} else {
2290	// more complicated
2291	colOrdered_ = ! colOrdered_;
2292	minorDim_ = rhs.majorDim_;
2293	majorDim_ = rhs.minorDim_;
2294	maxMajorDim_ = majorDim_ + extraForMajor;
2295	maxSize_ = CoinMax(size_ + extraElements,`1`);
2296	assert (maxMajorDim_>`0`);
2297	length_ = new int[maxMajorDim_];
2298	start_ = new CoinBigIndex[maxMajorDim_+`1`];
2299	element_ = new double[maxSize_];
2300	index_ = new int[maxSize_];
2301	bool hasGaps = rhs.size_<rhs.start_[rhs.majorDim_];
2302	CoinZeroN(length_, majorDim_);
2303	int i;
2304	if (hasGaps) {
2305	// has gaps
2306	for (i = `0`; i <rhs.majorDim_ ; ++i) {
2307	const CoinBigIndex first = rhs.start_[i];
2308	const CoinBigIndex last = first + rhs.length_[i];
2309	for (CoinBigIndex j = first; j < last; ++j) {
2310	assert( rhs.index_[j] < rhs.minorDim_ && rhs.index_[j]>=`0`);
2311	++length_[rhs.index_[j]];
2312	}
2313	}
2314	} else {
2315	// no gaps
2316	const CoinBigIndex last = rhs.start_[rhs.majorDim_];
2317	for (CoinBigIndex j = `0`; j < last; ++j) {
2318	assert( rhs.index_[j] < rhs.minorDim_ && rhs.index_[j]>=`0`);
2319	++length_[rhs.index_[j]];
2320	}
2321	}
2322	// Now do starts
2323	CoinBigIndex size=`0`;
2324	for (i = `0`; i <majorDim_ ; ++i) {
2325	start_[i]=size;
2326	size += length_[i];
2327	}
2328	start_[majorDim_]=size;
2329	assert (size==size_);
2330	for (i = `0`; i <rhs.majorDim_ ; ++i) {
2331	const CoinBigIndex first = rhs.start_[i];
2332	const CoinBigIndex last = first + rhs.length_[i];
2333	for (CoinBigIndex j = first; j < last; ++j) {
2334	const int ind = rhs.index_[j];
2335	CoinBigIndex put = start_[ind];
2336	start_[ind] = put +`1`;
2337	element_[put] = rhs.element_[j];
2338	index_[put] = i;
2339	}
2340	}
2341	// and re-adjust start_
2342	for (i = `0`; i < majorDim_; ++i) {
2343	start_[i] -= length_[i];
2344	}
2345	}
2346	}
2347	// Subset constructor (without gaps)
2348	CoinPackedMatrix::CoinPackedMatrix (const CoinPackedMatrix & rhs,
2349	int numberRows, const int * whichRow,
2350	int numberColumns,
2351	const int * whichColumn) :
2352	colOrdered_(true),
2353	extraGap_(`0.0`),
2354	extraMajor_(`0.0`),
2355	element_(NULL),
2356	index_(NULL),
2357	start_(NULL),
2358	length_(NULL),
2359	majorDim_(`0`),
2360	minorDim_(`0`),
2361	size_(`0`),
2362	maxMajorDim_(`0`),
2363	maxSize_(`0`)
2364	{
2365	if (numberRows<=`0`\|\|numberColumns<=`0`) {
2366	start_ = new CoinBigIndex[`1`];
2367	start_[`0`] = `0`;
2368	} else {
2369	if (!rhs.colOrdered_) {
2370	// just swap lists
2371	colOrdered_=false;
2372	const int * temp = whichRow;
2373	whichRow = whichColumn;
2374	whichColumn = temp;
2375	int n = numberRows;
2376	numberRows = numberColumns;
2377	numberColumns = n;
2378	}
2379	const double * element1 = rhs.element_;
2380	const int * index1 = rhs.index_;
2381	const CoinBigIndex * start1 = rhs.start_;
2382	const int * length1 = rhs.length_;
2383
2384	majorDim_ = numberColumns;
2385	maxMajorDim_ = numberColumns;
2386	minorDim_ = numberRows;
2387	// Throw exception if rhs empty
2388	if (rhs.majorDim_ <= `0` \|\| rhs.minorDim_ <= `0`)
2389	throw CoinError ("empty rhs", "subset constructor", "CoinPackedMatrix");
2390	// Array to say if an old row is in new copy
2391	int * newRow = new int [rhs.minorDim_];
2392	int iRow;
2393	for (iRow=`0`;iRow<rhs.minorDim_;iRow++)
2394	newRow[iRow] = -`1`;
2395	// and array for duplicating rows
2396	int * duplicateRow = new int [numberRows];
2397	int numberBad=`0`;
2398	int numberDuplicate=`0`;
2399	for (iRow=`0`;iRow<numberRows;iRow++) {
2400	duplicateRow[iRow] = -`1`;
2401	int kRow = whichRow[iRow];
2402	if (kRow>=`0` && kRow < rhs.minorDim_) {
2403	if (newRow[kRow]<`0`) {
2404	// first time
2405	newRow[kRow]=iRow;
2406	} else {
2407	// duplicate
2408	numberDuplicate++;
2409	int lastRow = newRow[kRow];
2410	newRow[kRow]=iRow;
2411	duplicateRow[iRow] = lastRow;
2412	}
2413	} else {
2414	// bad row
2415	numberBad++;
2416	}
2417	}
2418
2419	if (numberBad)
2420	throw CoinError ("bad minor entries",
2421	"subset constructor", "CoinPackedMatrix");
2422	// now get size and check columns
2423	size_ = `0`;
2424	int iColumn;
2425	numberBad=`0`;
2426	if (!numberDuplicate) {
2427	// No duplicates so can do faster
2428	// If not much smaller then use original size
2429	if (`3`majorDim_>`2`rhs.majorDim_&&
2430	`3`minorDim_>`2`rhs.minorDim_) {
2431	// now create arrays
2432	maxSize_=CoinMax(static_cast<CoinBigIndex> (`1`),rhs.size_);
2433	start_ = new CoinBigIndex [numberColumns+`1`];
2434	length_ = new int [numberColumns];
2435	index_ = new int[maxSize_];
2436	element_ = new double [maxSize_];
2437	// and fill them
2438	size_ = `0`;
2439	start_[`0`]=`0`;
2440	for (iColumn=`0`;iColumn<numberColumns;iColumn++) {
2441	int kColumn = whichColumn[iColumn];
2442	if (kColumn>=`0` && kColumn <rhs.majorDim_) {
2443	CoinBigIndex i;
2444	for (i=start1[kColumn];i<start1[kColumn]+length1[kColumn];i++) {
2445	int kRow = index1[i];
2446	double value = element1[i];
2447	kRow = newRow[kRow];
2448	if (kRow>=`0`) {
2449	index_[size_] = kRow;
2450	element_[size_++] = value;
2451	}
2452	}
2453	} else {
2454	// bad column
2455	numberBad++;
2456	}
2457	start_[iColumn+`1`] = size_;
2458	length_[iColumn] = size_ - start_[iColumn];
2459	}
2460	if (numberBad)
2461	throw CoinError ("bad major entries",
2462	"subset constructor", "CoinPackedMatrix");
2463	} else {
2464	for (iColumn=`0`;iColumn<numberColumns;iColumn++) {
2465	int kColumn = whichColumn[iColumn];
2466	if (kColumn>=`0` && kColumn <rhs.majorDim_) {
2467	CoinBigIndex i;
2468	for (i=start1[kColumn];i<start1[kColumn]+length1[kColumn];i++) {
2469	int kRow = index1[i];
2470	kRow = newRow[kRow];
2471	if (kRow>=`0`)
2472	size_++;
2473	}
2474	} else {
2475	// bad column
2476	numberBad++;
2477	}
2478	}
2479	if (numberBad)
2480	throw CoinError ("bad major entries",
2481	"subset constructor", "CoinPackedMatrix");
2482	// now create arrays
2483	maxSize_=CoinMax(static_cast<CoinBigIndex> (`1`),size_);
2484	start_ = new CoinBigIndex [numberColumns+`1`];
2485	length_ = new int [numberColumns];
2486	index_ = new int[maxSize_];
2487	element_ = new double [maxSize_];
2488	// and fill them
2489	size_ = `0`;
2490	start_[`0`]=`0`;
2491	for (iColumn=`0`;iColumn<numberColumns;iColumn++) {
2492	int kColumn = whichColumn[iColumn];
2493	CoinBigIndex i;
2494	for (i=start1[kColumn];i<start1[kColumn]+length1[kColumn];i++) {
2495	int kRow = index1[i];
2496	double value = element1[i];
2497	kRow = newRow[kRow];
2498	if (kRow>=`0`) {
2499	index_[size_] = kRow;
2500	element_[size_++] = value;
2501	}
2502	}
2503	start_[iColumn+`1`] = size_;
2504	length_[iColumn] = size_ - start_[iColumn];
2505	}
2506	}
2507	} else {
2508	for (iColumn=`0`;iColumn<numberColumns;iColumn++) {
2509	int kColumn = whichColumn[iColumn];
2510	if (kColumn>=`0` && kColumn <rhs.majorDim_) {
2511	CoinBigIndex i;
2512	for (i=start1[kColumn];i<start1[kColumn]+length1[kColumn];i++) {
2513	int kRow = index1[i];
2514	kRow = newRow[kRow];
2515	while (kRow>=`0`) {
2516	size_++;
2517	kRow = duplicateRow[kRow];
2518	}
2519	}
2520	} else {
2521	// bad column
2522	numberBad++;
2523	}
2524	}
2525	if (numberBad)
2526	throw CoinError ("bad major entries",
2527	"subset constructor", "CoinPackedMatrix");
2528	// now create arrays
2529	maxSize_=CoinMax(static_cast<CoinBigIndex> (`1`),size_);
2530	start_ = new CoinBigIndex [numberColumns+`1`];
2531	length_ = new int [numberColumns];
2532	index_ = new int[maxSize_];
2533	element_ = new double [maxSize_];
2534	// and fill them
2535	size_ = `0`;
2536	start_[`0`]=`0`;
2537	for (iColumn=`0`;iColumn<numberColumns;iColumn++) {
2538	int kColumn = whichColumn[iColumn];
2539	CoinBigIndex i;
2540	for (i=start1[kColumn];i<start1[kColumn]+length1[kColumn];i++) {
2541	int kRow = index1[i];
2542	double value = element1[i];
2543	kRow = newRow[kRow];
2544	while (kRow>=`0`) {
2545	index_[size_] = kRow;
2546	element_[size_++] = value;
2547	kRow = duplicateRow[kRow];
2548	}
2549	}
2550	start_[iColumn+`1`] = size_;
2551	length_[iColumn] = size_ - start_[iColumn];
2552	}
2553	}
2554	delete [] newRow;
2555	delete [] duplicateRow;
2556	}
2557	}
2558
2559
2560	//-----------------------------------------------------------------------------
2561
2562	CoinPackedMatrix::~CoinPackedMatrix ()
2563	{
2564	gutsOfDestructor();
2565	}
2566
2567	//#############################################################################
2568	//#############################################################################
2569	//#############################################################################
2570
2571	void
2572	CoinPackedMatrix::gutsOfDestructor()
2573	{
2574	delete[] length_;
2575	delete[] start_;
2576	delete[] index_;
2577	delete[] element_;
2578	length_ = `0`;
2579	start_ = `0`;
2580	index_ = `0`;
2581	element_ = `0`;
2582	}
2583
2584	//#############################################################################
2585
2586	void
2587	CoinPackedMatrix::gutsOfCopyOf(const bool colordered,
2588	const int minor, const int major,
2589	const CoinBigIndex numels,
2590	const double * elem, const int * ind,
2591	const CoinBigIndex * start, const int * len,
2592	const double extraMajor, const double extraGap)
2593	{
2594	colOrdered_ = colordered;
2595	majorDim_ = major;
2596	minorDim_ = minor;
2597	size_ = numels;
2598
2599	extraGap_ = extraGap;
2600	extraMajor_ = extraMajor;
2601
2602	maxMajorDim_ = CoinLengthWithExtra(majorDim_, extraMajor_);
2603
2604	if (maxMajorDim_ > `0`) {
2605	delete [] length_;
2606	length_ = new int[maxMajorDim_];
2607	if (len == `0`) {
2608	std::adjacent_difference(start + `1`, start + (major + `1`), length_);
2609	length_[`0`] -= start[`0`];
2610	} else {
2611	CoinMemcpyN(len, major, length_);
2612	}
2613	delete [] start_;
2614	start_ = new CoinBigIndex[maxMajorDim_+`1`];
2615	start_[`0`]=`0`;
2616	CoinMemcpyN(start, major+`1`, start_);
2617	} else {
2618	// empty but be safe
2619	delete [] length_;
2620	length_ = NULL;
2621	delete [] start_;
2622	start_ = new CoinBigIndex[`1`];
2623	start_[`0`]=`0`;
2624	}
2625
2626	maxSize_ = maxMajorDim_ > `0` ? start_[major] : `0`;
2627	maxSize_ = CoinLengthWithExtra(maxSize_, extraMajor_);
2628
2629	if (maxSize_ > `0`) {
2630	delete [] element_;
2631	delete []index_;
2632	element_ = new double[maxSize_];
2633	index_ = new int[maxSize_];
2634	// we can't just simply memcpy these content over, because that can
2635	// upset memory debuggers like purify if there were gaps and those gaps
2636	// were uninitialized memory blocks
2637	for (int i = majorDim_ - `1`; i >= `0`; --i) {
2638	CoinMemcpyN(ind + start[i], length_[i], index_ + start_[i]);
2639	CoinMemcpyN(elem + start[i], length_[i], element_ + start_[i]);
2640	}
2641	}
2642	}
2643
2644	//#############################################################################
2645
2646	void
2647	CoinPackedMatrix::gutsOfCopyOfNoGaps(const bool colordered,
2648	const int minor, const int major,
2649	const double * elem, const int * ind,
2650	const CoinBigIndex * start)
2651	{
2652	colOrdered_ = colordered;
2653	majorDim_ = major;
2654	minorDim_ = minor;
2655	size_ = start[majorDim_];
2656
2657	extraGap_ = `0`;
2658	extraMajor_ = `0`;
2659
2660	maxMajorDim_ = majorDim_;
2661
2662	// delete all arrays
2663	delete [] length_;
2664	delete [] start_;
2665	delete [] element_;
2666	delete [] index_;
2667
2668	if (maxMajorDim_ > `0`) {
2669	length_ = new int[maxMajorDim_];
2670	assert (!start[`0`]);
2671	start_ = new CoinBigIndex[maxMajorDim_+`1`];
2672	start_[`0`]=`0`;
2673	CoinBigIndex last = `0`;
2674	for (int i=`0`;i<majorDim_;i++) {
2675	CoinBigIndex first = last;
2676	last = start[i+`1`];
2677	length_[i] = last-first;
2678	start_[i+`1`]=last;
2679	}
2680	} else {
2681	// empty but be safe
2682	length_ = NULL;
2683	start_ = new CoinBigIndex[`1`];
2684	start_[`0`]=`0`;
2685	}
2686
2687	maxSize_ = start_[majorDim_];
2688
2689	if (maxSize_ > `0`) {
2690	element_ = new double[maxSize_];
2691	index_ = new int[maxSize_];
2692	CoinMemcpyN(ind , maxSize_, index_);
2693	CoinMemcpyN(elem , maxSize_, element_);
2694	} else {
2695	element_ = NULL;
2696	index_ = NULL;
2697	}
2698	}
2699
2700	//#############################################################################
2701
2702	void
2703	CoinPackedMatrix::gutsOfOpEqual(const bool colordered,
2704	const int minor, const int major,
2705	const CoinBigIndex numels,
2706	const double * elem, const int * ind,
2707	const CoinBigIndex * start, const int * len)
2708	{
2709	colOrdered_ = colordered;
2710	majorDim_ = major;
2711	minorDim_ = minor;
2712	size_ = numels;
2713	if (!len && numels > `0` && numels==start[major] && start[`0`]==`0`) {
2714	// No gaps - do faster
2715	if (major>maxMajorDim_\|\|!start_) {
2716	maxMajorDim_ = major;
2717	delete [] length_;
2718	length_ = new int[maxMajorDim_];
2719	delete [] start_;
2720	start_ = new CoinBigIndex[maxMajorDim_+`1`];
2721	}
2722	CoinMemcpyN(start,major+`1`,start_);
2723	std::adjacent_difference(start + `1`, start + (major + `1`), length_);
2724	if (numels>maxSize_\|\|!element_) {
2725	maxSize_=numels;
2726	delete [] element_;
2727	delete [] index_;
2728	element_ = new double[maxSize_];
2729	index_ = new int[maxSize_];
2730	}
2731	CoinMemcpyN(ind,numels,index_);
2732	CoinMemcpyN(elem,numels,element_);
2733	} else {
2734
2735	maxMajorDim_ = CoinLengthWithExtra(majorDim_, extraMajor_);
2736
2737	int i;
2738	if (maxMajorDim_ > `0`) {
2739	delete [] length_;
2740	length_ = new int[maxMajorDim_];
2741	if (len == `0`) {
2742	std::adjacent_difference(start + `1`, start + (major + `1`), length_);
2743	length_[`0`] -= start[`0`];
2744	} else {
2745	CoinMemcpyN(len, major, length_);
2746	}
2747	delete [] start_;
2748	start_ = new CoinBigIndex[maxMajorDim_+`1`];
2749	start_[`0`] = `0`;
2750	if (extraGap_ == `0`) {
2751	for (i = `0`; i < major; ++i)
2752	start_[i+`1`] = start_[i] + length_[i];
2753	} else {
2754	const double extra_gap = extraGap_;
2755	for (i = `0`; i < major; ++i)
2756	start_[i+`1`] = start_[i] + CoinLengthWithExtra(length_[i], extra_gap);
2757	}
2758	} else {
2759	// empty matrix
2760	delete [] start_;
2761	start_ = new CoinBigIndex[`1`];
2762	start_[`0`] = `0`;
2763	}
2764
2765	maxSize_ = maxMajorDim_ > `0` ? start_[major] : `0`;
2766	maxSize_ = CoinLengthWithExtra(maxSize_, extraMajor_);
2767
2768	if (maxSize_ > `0`) {
2769	delete [] element_;
2770	delete [] index_;
2771	element_ = new double[maxSize_];
2772	index_ = new int[maxSize_];
2773	assert (maxSize_>=start_[majorDim_-`1`]+length_[majorDim_-`1`]);
2774	// we can't just simply memcpy these content over, because that can
2775	// upset memory debuggers like purify if there were gaps and those gaps
2776	// were uninitialized memory blocks
2777	for (i = majorDim_ - `1`; i >= `0`; --i) {
2778	CoinMemcpyN(ind + start[i], length_[i], index_ + start_[i]);
2779	CoinMemcpyN(elem + start[i], length_[i], element_ + start_[i]);
2780	}
2781	}
2782	}
2783	#ifndef NDEBUG
2784	for (int i = majorDim_ - `1`; i >= `0`; --i) {
2785	const CoinBigIndex last = getVectorLast(i);
2786	for (CoinBigIndex j = getVectorFirst(i); j < last; ++j) {
2787	int index = index_[j];
2788	assert (index>=`0`&&index<minorDim_);
2789	}
2790	}
2791	#endif
2792	}
2793
2794	//#############################################################################
2795
2796	// This routine is called only if we MUST resize!
2797	void
2798	CoinPackedMatrix::resizeForAddingMajorVectors(const int numVec,
2799	const int * lengthVec)
2800	{
2801	const double extra_gap = extraGap_;
2802	int i;
2803
2804	maxMajorDim_ =
2805	CoinMax(maxMajorDim_, CoinLengthWithExtra(majorDim_ + numVec, extraMajor_));
2806
2807	CoinBigIndex * newStart = new CoinBigIndex[maxMajorDim_ + `1`];
2808	int * newLength = new int[maxMajorDim_];
2809
2810	CoinMemcpyN(length_, majorDim_, newLength);
2811	// fake that the new vectors are there
2812	CoinMemcpyN(lengthVec, numVec, newLength + majorDim_);
2813	majorDim_ += numVec;
2814
2815	newStart[`0`] = `0`;
2816	if (extra_gap == `0`) {
2817	for (i = `0`; i < majorDim_; ++i)
2818	newStart[i+`1`] = newStart[i] + newLength[i];
2819	} else {
2820	for (i = `0`; i < majorDim_; ++i)
2821	newStart[i+`1`] = newStart[i] + CoinLengthWithExtra(newLength[i],extra_gap);
2822	}
2823
2824	maxSize_ =
2825	CoinMax(maxSize_, CoinLengthWithExtra(newStart[majorDim_], extraMajor_));
2826	majorDim_ -= numVec;
2827
2828	int * newIndex = new int[maxSize_];
2829	double * newElem = new double[maxSize_];
2830	for (i = majorDim_ - `1`; i >= `0`; --i) {
2831	CoinMemcpyN(index_ + start_[i], length_[i], newIndex + newStart[i]);
2832	CoinMemcpyN(element_ + start_[i], length_[i], newElem + newStart[i]);
2833	}
2834
2835	gutsOfDestructor();
2836	start_ = newStart;
2837	length_ = newLength;
2838	index_ = newIndex;
2839	element_ = newElem;
2840	}
2841
2842
2843	//#############################################################################
2844
2845	void
2846	CoinPackedMatrix::resizeForAddingMinorVectors(const int * addedEntries)
2847	{
2848	int i;
2849	maxMajorDim_ =
2850	CoinMax(CoinLengthWithExtra(majorDim_, extraMajor_), maxMajorDim_);
2851	CoinBigIndex * newStart = new CoinBigIndex[maxMajorDim_ + `1`];
2852	int * newLength = new int[maxMajorDim_];
2853	// increase the lengths temporarily so that the correct new start positions
2854	// can be easily computed (it's faster to modify the lengths and reset them
2855	// than do a test for every entry when the start positions are computed.
2856	for (i = majorDim_ - `1`; i >= `0`; --i)
2857	newLength[i] = length_[i] + addedEntries[i];
2858
2859	newStart[`0`] = `0`;
2860	if (extraGap_ == `0`) {
2861	for (i = `0`; i < majorDim_; ++i)
2862	newStart[i+`1`] = newStart[i] + newLength[i];
2863	} else {
2864	const double eg = extraGap_;
2865	for (i = `0`; i < majorDim_; ++i)
2866	newStart[i+`1`] = newStart[i] + CoinLengthWithExtra(newLength[i], eg);
2867	}
2868
2869	// reset the lengths
2870	for (i = majorDim_ - `1`; i >= `0`; --i)
2871	newLength[i] -= addedEntries[i];
2872
2873	maxSize_ =
2874	CoinMax(maxSize_, CoinLengthWithExtra(newStart[majorDim_], extraMajor_));
2875	int * newIndex = new int[maxSize_];
2876	double * newElem = new double[maxSize_];
2877	for (i = majorDim_ - `1`; i >= `0`; --i) {
2878	CoinMemcpyN(index_ + start_[i], length_[i],
2879	newIndex + newStart[i]);
2880	CoinMemcpyN(element_ + start_[i], length_[i],
2881	newElem + newStart[i]);
2882	}
2883
2884	gutsOfDestructor();
2885	start_ = newStart;
2886	length_ = newLength;
2887	index_ = newIndex;
2888	element_ = newElem;
2889	}
2890
2891	//#############################################################################
2892	//#############################################################################
2893
2894	void
2895	CoinPackedMatrix::dumpMatrix(const char* fname) const
2896	{
2897	if (! fname) {
2898	printf("Dumping matrix...\n\n");
2899	printf("colordered: %i\n", isColOrdered() ? `1` : `0`);
2900	const int major = getMajorDim();
2901	const int minor = getMinorDim();
2902	printf("major: %i minor: %i\n", major, minor);
2903	for (int i = `0`; i < major; ++i) {
2904	printf("vec %i has length %i with entries:\n", i, length_[i]);
2905	for (CoinBigIndex j = start_[i]; j < start_[i] + length_[i]; ++j) {
2906	printf(" %15i %40.25f\n", index_[j], element_[j]);
2907	}
2908	}
2909	printf("\nFinished dumping matrix\n");
2910	} else {
2911	FILE* out = fopen(fname, "w");
2912	fprintf(out, "Dumping matrix...\n\n");
2913	fprintf(out, "colordered: %i\n", isColOrdered() ? `1` : `0`);
2914	const int major = getMajorDim();
2915	const int minor = getMinorDim();
2916	fprintf(out, "major: %i minor: %i\n", major, minor);
2917	for (int i = `0`; i < major; ++i) {
2918	fprintf(out, "vec %i has length %i with entries:\n", i, length_[i]);
2919	for (CoinBigIndex j = start_[i]; j < start_[i] + length_[i]; ++j) {
2920	fprintf(out, " %15i %40.25f\n", index_[j], element_[j]);
2921	}
2922	}
2923	fprintf(out, "\nFinished dumping matrix\n");
2924	fclose(out);
2925	}
2926	}
2927	void
2928	CoinPackedMatrix::printMatrixElement (const int row_val,
2929	const int col_val) const
2930	{
2931	int major_index, minor_index;
2932	if (isColOrdered()) {
2933	major_index = col_val;
2934	minor_index = row_val;
2935	} else {
2936	major_index = row_val;
2937	minor_index = col_val;
2938	}
2939	if (major_index < `0` \|\| major_index > getMajorDim()-`1`) {
2940	std::cout
2941	<< "Major index " << major_index << " not in range 0.."
2942	<< getMajorDim()-`1` << std::endl ;
2943	} else if (minor_index < `0` \|\| minor_index > getMinorDim()-`1`) {
2944	std::cout
2945	<< "Minor index " << minor_index << " not in range 0.."
2946	<< getMinorDim()-`1` << std::endl ;
2947	} else {
2948	CoinBigIndex curr_point = start_[major_index];
2949	const CoinBigIndex stop_point = curr_point+length_[major_index];
2950	double aij = `0.0` ;
2951	for ( ; curr_point < stop_point ; curr_point++) {
2952	if (index_[curr_point] == minor_index) {
2953	aij = element_[curr_point];
2954	break;
2955	}
2956	}
2957	std::cout << aij ;
2958	}
2959	}
2960	#ifndef CLP_NO_VECTOR
2961	bool
2962	CoinPackedMatrix::isEquivalent2(const CoinPackedMatrix& rhs) const
2963	{
2964	CoinRelFltEq eq;
2965	// Both must be column order or both row ordered and must be of same size
2966	if (isColOrdered() ^ rhs.isColOrdered()) {
2967	std::cerr <<"Ordering "<<isColOrdered()<<
2968	" rhs - "<<rhs.isColOrdered()<<std::endl;
2969	return false;
2970	}
2971	if (getNumCols() != rhs.getNumCols()) {
2972	std::cerr <<"NumCols "<<getNumCols()<<
2973	" rhs - "<<rhs.getNumCols()<<std::endl;
2974	return false;
2975	}
2976	if (getNumRows() != rhs.getNumRows()) {
2977	std::cerr <<"NumRows "<<getNumRows()<<
2978	" rhs - "<<rhs.getNumRows()<<std::endl;
2979	return false;
2980	}
2981	if (getNumElements() != rhs.getNumElements()) {
2982	std::cerr <<"NumElements "<<getNumElements()<<
2983	" rhs - "<<rhs.getNumElements()<<std::endl;
2984	return false;
2985	}
2986
2987	for (int i=getMajorDim()-`1`; i >= `0`; --i) {
2988	CoinShallowPackedVector pv = getVector(i);
2989	CoinShallowPackedVector rhsPv = rhs.getVector(i);
2990	if ( !pv.isEquivalent(rhsPv,eq) ) {
2991	std::cerr <<"vector # "<<i <<" nel "<<pv.getNumElements()<<
2992	" rhs - "<<rhsPv.getNumElements()<<std::endl;
2993	int j;
2994	const int * inds = pv.getIndices();
2995	const double * elems = pv.getElements();
2996	const int * inds2 = rhsPv.getIndices();
2997	const double * elems2 = rhsPv.getElements();
2998	for ( j = `0` ;j < pv.getNumElements() ; ++j) {
2999	double diff = elems[j]-elems2[j];
3000	if (diff) {
3001	std::cerr <<j <<"( "<<inds[j]<<", "<<elems[j]<<"), rhs ( "<<
3002	inds2[j]<<", "<<elems2[j]<<") diff "<<
3003	diff <<std::endl;
3004	const int * xx = reinterpret_cast<const int *> (elems+j);
3005	printf("%x %x",xx[`0`],xx[`1`]);
3006	xx = reinterpret_cast<const int *> (elems2+j);
3007	printf(" %x %x\n",xx[`0`],xx[`1`]);
3008	}
3009	}
3010	//return false;
3011	}
3012	}
3013	return true;
3014	}
3015	#else
3016	/ Equivalence.*
3017	Two matrices are equivalent if they are both by rows or both by columns,
3018	they have the same dimensions, and each vector is equivalent.
3019	In this method the FloatEqual function operator can be specified.
3020	*/
3021	bool
3022	CoinPackedMatrix::isEquivalent(const CoinPackedMatrix& rhs, const CoinRelFltEq& eq) const
3023	{
3024	// Both must be column order or both row ordered and must be of same size
3025	if ((isColOrdered() ^ rhs.isColOrdered()) \|\|
3026	(getNumCols() != rhs.getNumCols()) \|\|
3027	(getNumRows() != rhs.getNumRows()) \|\|
3028	(getNumElements() != rhs.getNumElements()))
3029	return false;
3030
3031	const int major = getMajorDim();
3032	const int minor = getMinorDim();
3033	double * values = new double[minor];
3034	memset(values,`0`,minor*sizeof(double));
3035	bool same=true;
3036	for (int i = `0`; i < major; ++i) {
3037	int length = length_[i];
3038	if (length!=rhs.length_[i]) {
3039	same=false;
3040	break;
3041	} else {
3042	CoinBigIndex j;
3043	for ( j = start_[i]; j < start_[i] + length; ++j) {
3044	int index = index_[j];
3045	values[index]=element_[j];
3046	}
3047	for ( j = rhs.start_[i]; j < rhs.start_[i] + length; ++j) {
3048	int index = index_[j];
3049	double oldValue = values[index];
3050	values[index]=`0.0`;
3051	if (!eq(oldValue,rhs.element_[j])) {
3052	same=false;
3053	break;
3054	}
3055	}
3056	if (!same)
3057	break;
3058	}
3059	}
3060	delete [] values;
3061	return same;
3062	}
3063	#endif
3064	bool CoinPackedMatrix::isEquivalent(const CoinPackedMatrix& rhs) const
3065	{
3066	return isEquivalent(rhs,CoinRelFltEq ());
3067	}
3068	/ Sort all columns so indices are increasing.in each column /
3069	void
3070	CoinPackedMatrix::orderMatrix()
3071	{
3072	for (int i=`0`;i<majorDim_;i++) {
3073	CoinBigIndex start = start_[i];
3074	CoinBigIndex end = start + length_[i];
3075	CoinSort_2(index_+start,index_+end,element_+start);
3076	}
3077	}
3078	/ Append a set of rows/columns to the end of the matrix. Returns number of errors*
3079	i.e. if any of the new rows/columns contain an index that's larger than the
3080	number of columns-1/rows-1 (if numberOther>0) or duplicates
3081	This version is easy one i.e. adding columns to column ordered /*
3082	int
3083	CoinPackedMatrix::appendMajor(const int number,
3084	const CoinBigIndex * starts, const int * index,
3085	const double * element, int numberOther)
3086	{
3087	int i;
3088	int numberErrors=`0`;
3089	CoinBigIndex numberElements = starts[number];
3090	if (majorDim_ + number > maxMajorDim_ \|\|
3091	getLastStart() + numberElements > maxSize_) {
3092	// we got to resize before we add. note that the resizing method
3093	// properly fills out start_ and length_ for the major-dimension
3094	// vectors to be added!
3095	if (!extraGap_&&!extraMajor_&&numberOther<=`0`&&!hasGaps()) {
3096	// can do faster
3097	if (majorDim_+number>maxMajorDim_) {
3098	maxMajorDim_ = majorDim_+number;
3099	int * newLength = new int[maxMajorDim_];
3100	CoinMemcpyN(length_, majorDim_, newLength);
3101	delete [] length_;
3102	length_ = newLength;
3103	CoinBigIndex * newStart = new CoinBigIndex[maxMajorDim_ + `1`];
3104	CoinMemcpyN(start_, majorDim_+`1`, newStart);
3105	delete [] start_;
3106	start_ = newStart;
3107	}
3108	if (size_+numberElements>maxSize_) {
3109	maxSize_ = size_+numberElements;
3110	double * newElem = new double[maxSize_];
3111	CoinMemcpyN(element_,size_,newElem);
3112	delete [] element_;
3113	element_ = newElem;
3114	int * newIndex = new int[maxSize_];
3115	CoinMemcpyN(index_,size_,newIndex);
3116	delete [] index_;
3117	index_ = newIndex;
3118	}
3119	CoinMemcpyN(index,numberElements,index_+size_);
3120	CoinMemcpyN(element,numberElements,element_+size_);
3121	i=majorDim_;
3122	starts -= majorDim_;
3123	majorDim_ += number;
3124	int iStart=`0`;
3125	for (;i<majorDim_;i++) {
3126	int next = starts[i+`1`];
3127	int length = next-iStart;
3128	length_[i]=length;
3129	iStart=next;
3130	size_ += length;
3131	start_[i+`1`]=size_;
3132	}
3133	return `0`;
3134	} else {
3135	int * length = new int[number];
3136	for (i=`0`;i<number;i++)
3137	length[i]=starts[i+`1`]-starts[i];
3138	resizeForAddingMajorVectors(number, length);
3139	delete [] length;
3140	}
3141	if (numberOther>`0`) {
3142	char * which = new char[numberOther];
3143	memset(which,`0`,numberOther);
3144	for (i = `0`; i < number; i++) {
3145	CoinBigIndex put = start_[majorDim_+i];
3146	CoinBigIndex j;
3147	for ( j=starts[i];j<starts[i+`1`];j++) {
3148	int iIndex = index[j];
3149	element_[put]=element[j];
3150	if (iIndex>=`0`&&iIndex<numberOther) {
3151	if (!which[iIndex])
3152	which[iIndex]=`1`;
3153	else
3154	numberErrors++;
3155	} else {
3156	numberErrors++;
3157	}
3158	index_[put++]=iIndex;
3159	}
3160	for ( j=starts[i];j<starts[i+`1`];j++) {
3161	int iIndex = index[j];
3162	if (iIndex>=`0`&&iIndex<numberOther)
3163	which[iIndex]=`0`;
3164	}
3165	}
3166	delete [] which;
3167	} else {
3168	// easy
3169	int lastMinor=-`1`;
3170	if (!extraGap_) {
3171	// just one copy
3172	int * index2 = index_+start_[majorDim_];
3173	for (CoinBigIndex j=`0`;j<numberElements;j++) {
3174	int iIndex = index[j];
3175	index2[j] = iIndex;
3176	lastMinor = CoinMax(lastMinor,iIndex);
3177	}
3178	CoinMemcpyN(element,numberElements,element_+start_[majorDim_]);
3179	} else {
3180	start_ += majorDim_;
3181	for (i = `0`; i < number; i++) {
3182	int length = starts[i+`1`]-starts[i];
3183	int * index2 = index_+start_[i];
3184	const int * index1 = index+starts[i];
3185	for (CoinBigIndex j=`0`;j<length;j++) {
3186	int iIndex = index1[j];
3187	index2[j] = iIndex;
3188	lastMinor = CoinMax(lastMinor,iIndex);
3189	}
3190	CoinMemcpyN(element + starts[i], length,
3191	element_ + start_[i]);
3192	}
3193	start_ -= majorDim_;
3194	}
3195	// update minorDim if necessary
3196	minorDim_ = CoinMax(minorDim_,lastMinor+`1`);
3197	}
3198	} else {
3199	if (numberOther>`0`) {
3200	char * which = new char[numberOther];
3201	memset(which,`0`,numberOther);
3202	for (i = `0`; i < number; i++) {
3203	CoinBigIndex put = start_[majorDim_+i];
3204	CoinBigIndex j;
3205	for ( j=starts[i];j<starts[i+`1`];j++) {
3206	int iIndex = index[j];
3207	element_[put]=element[j];
3208	if (iIndex>=`0`&&iIndex<numberOther) {
3209	if (!which[iIndex])
3210	which[iIndex]=`1`;
3211	else
3212	numberErrors++;
3213	} else {
3214	numberErrors++;
3215	}
3216	index_[put++]=iIndex;
3217	}
3218	start_[majorDim_+i+`1`] = put;
3219	length_[majorDim_+i] = put-start_[majorDim_+i];
3220	for ( j=starts[i];j<starts[i+`1`];j++) {
3221	int iIndex = index[j];
3222	if (iIndex>=`0`&&iIndex<numberOther)
3223	which[iIndex]=`0`;
3224	}
3225	}
3226	delete [] which;
3227	} else {
3228	// easy
3229	int lastMinor=-`1`;
3230	if (!extraGap_) {
3231	// just one copy
3232	// just one copy
3233	int * index2 = index_+start_[majorDim_];
3234	for (CoinBigIndex j=`0`;j<numberElements;j++) {
3235	int iIndex = index[j];
3236	index2[j] = iIndex;
3237	lastMinor = CoinMax(lastMinor,iIndex);
3238	}
3239	CoinMemcpyN(element,numberElements,element_+start_[majorDim_]);
3240	start_ += majorDim_;
3241	for (i = `0`; i < number; i++) {
3242	int length = starts[i+`1`]-starts[i];
3243	start_[i+`1`] = start_[i] + length;
3244	length_[majorDim_+i] = length;
3245	}
3246	start_ -= majorDim_;
3247	} else {
3248	start_ += majorDim_;
3249	for (i = `0`; i < number; i++) {
3250	int length = starts[i+`1`]-starts[i];
3251	int * index2 = index_+start_[i];
3252	const int * index1 = index+starts[i];
3253	for (CoinBigIndex j=`0`;j<length;j++) {
3254	int iIndex = index1[j];
3255	index2[j] = iIndex;
3256	lastMinor = CoinMax(lastMinor,iIndex);
3257	}
3258	CoinMemcpyN(element + starts[i], length,
3259	element_ + start_[i]);
3260	start_[i+`1`] = start_[i] + length;
3261	length_[majorDim_+i] = length;
3262	}
3263	start_ -= majorDim_;
3264	}
3265	// update minorDim if necessary
3266	minorDim_ = CoinMax(minorDim_,lastMinor+`1`);
3267	}
3268	}
3269	majorDim_ += number;
3270	size_ += numberElements;
3271	#ifndef NDEBUG
3272	int checkSize=`0`;
3273	for (int i=`0`;i<majorDim_;i++) {
3274	checkSize += length_[i];
3275	}
3276	assert (checkSize==size_);
3277	#endif
3278	return numberErrors;
3279	}
3280	/ Append a set of rows/columns to the end of the matrix. Returns number of errors*
3281	i.e. if any of the new rows/columns contain an index that's larger than the
3282	number of columns-1/rows-1 (if numberOther>0) or duplicates
3283	This version is harder one i.e. adding columns to row ordered /*
3284	int
3285	CoinPackedMatrix::appendMinor(const int number,
3286	const CoinBigIndex * starts, const int * index,
3287	const double * element, int numberOther)
3288	{
3289	int i;
3290	int numberErrors=`0`;
3291	// first compute how many entries will be added to each major-dimension
3292	// vector, and if needed, resize the matrix to accommodate all
3293	int * addedEntries = NULL;
3294	if (numberOther>`0`) {
3295	addedEntries = new int[majorDim_];
3296	CoinZeroN(addedEntries,majorDim_);
3297	numberOther=majorDim_;
3298	char * which = new char[numberOther];
3299	memset(which,`0`,numberOther);
3300	for (i = `0`; i < number; i++) {
3301	CoinBigIndex j;
3302	for ( j=starts[i];j<starts[i+`1`];j++) {
3303	int iIndex = index[j];
3304	if (iIndex>=`0`&&iIndex<numberOther) {
3305	addedEntries[iIndex]++;
3306	if (!which[iIndex])
3307	which[iIndex]=`1`;
3308	else
3309	numberErrors++;
3310	} else {
3311	numberErrors++;
3312	}
3313	}
3314	for ( j=starts[i];j<starts[i+`1`];j++) {
3315	int iIndex = index[j];
3316	if (iIndex>=`0`&&iIndex<numberOther)
3317	which[iIndex]=`0`;
3318	}
3319	}
3320	delete [] which;
3321	} else {
3322	int largest = majorDim_-`1`;
3323	for (i = `0`; i < number; i++) {
3324	CoinBigIndex j;
3325	for ( j=starts[i];j<starts[i+`1`];j++) {
3326	int iIndex = index[j];
3327	largest = CoinMax(largest,iIndex);
3328	}
3329	}
3330	if (largest+`1`>majorDim_) {
3331	if (isColOrdered())
3332	setDimensions(-`1`,largest+`1`);
3333	else
3334	setDimensions(largest+`1`,-`1`);
3335	}
3336	addedEntries = new int[majorDim_];
3337	CoinZeroN(addedEntries,majorDim_);
3338	// no checking
3339	for (i = `0`; i < number; i++) {
3340	CoinBigIndex j;
3341	for ( j=starts[i];j<starts[i+`1`];j++) {
3342	int iIndex = index[j];
3343	addedEntries[iIndex]++;
3344	}
3345	}
3346	}
3347	for (i = majorDim_ - `1`; i >= `0`; i--) {
3348	if (start_[i] + length_[i] + addedEntries[i] > start_[i+`1`])
3349	break;
3350	}
3351	if (i >= `0`)
3352	resizeForAddingMinorVectors(addedEntries);
3353	delete[] addedEntries;
3354
3355	// now insert the entries of matrix
3356	for (i = `0`; i < number; i++) {
3357	CoinBigIndex j;
3358	for ( j=starts[i];j<starts[i+`1`];j++) {
3359	int iIndex = index[j];
3360	element_[start_[iIndex] + length_[iIndex]] = element[j];
3361	index_[start_[iIndex] + (length_[iIndex]++)] = minorDim_;
3362	}
3363	++minorDim_;
3364	}
3365	size_ += starts[number];
3366	#ifndef NDEBUG
3367	int checkSize=`0`;
3368	for (int i=`0`;i<majorDim_;i++) {
3369	checkSize += length_[i];
3370	}
3371	assert (checkSize==size_);
3372	#endif
3373	return numberErrors;
3374	}
3375	//#define ADD_ROW_ANALYZE
3376	#ifdef ADD_ROW_ANALYZE
3377	static int xxxxxx[`10`]={`0`,`0`,`0`,`0`,`0`,`0`,`0`,`0`,`0`,`0`};
3378	#endif
3379	/ Append a set of rows/columns to the end of the matrix. This case is*
3380	when we know there are no gaps and majorDim_ will not change
3381	This version is harder one i.e. adding columns to row ordered /*
3382	void
3383	CoinPackedMatrix::appendMinorFast(const int number,
3384	const CoinBigIndex * starts, const int * index,
3385	const double * element)
3386	{
3387	#ifdef ADD_ROW_ANALYZE
3388	xxxxxx[`0`]++;
3389	#endif
3390	// first compute how many entries will be added to each major-dimension
3391	// vector, and if needed, resize the matrix to accommodate all
3392	// Will be used as new start array
3393	CoinBigIndex * newStart = new CoinBigIndex [maxMajorDim_+`1`];
3394	CoinZeroN(newStart,maxMajorDim_);
3395	// no checking
3396	int numberAdded = starts[number];
3397	for (CoinBigIndex j = `0`; j < numberAdded; j++) {
3398	int iIndex = index[j];
3399	newStart[iIndex]++;
3400	}
3401	int packType=`0`;
3402	#ifdef ADD_ROW_ANALYZE
3403	int nBad=`0`;
3404	#endif
3405	if (size_+numberAdded<=maxSize_) {
3406	CoinBigIndex nextStart=start_[majorDim_];
3407	// could do other way and then stop moving
3408	for (int i = majorDim_ - `1`; i >= `0`; i--) {
3409	CoinBigIndex start = start_[i];
3410	if (start + length_[i] + newStart[i] <= nextStart) {
3411	nextStart=start;
3412	} else {
3413	packType=-`1`;
3414	#ifdef ADD_ROW_ANALYZE
3415	nBad++;
3416	#else
3417	break;
3418	#endif
3419	}
3420	}
3421	} else {
3422	// Need more space
3423	packType=`1`;
3424	}
3425	#ifdef ADD_ROW_ANALYZE
3426	if (!hasGaps())
3427	xxxxxx[`9`]++;
3428	if (packType==-`1`&&nBad<`6`)
3429	packType=nBad+`1`;
3430	xxxxxx[packType+`2`]++;
3431	if ((xxxxxx[`0`]%`100`)==`0`) {
3432	printf("Append ");
3433	for (int i=`0`;i<`10`;i++)
3434	printf("%d ",xxxxxx[i]);
3435	printf("\n");
3436	}
3437	#endif
3438	if (hasGaps()&&packType)
3439	packType=`1`;
3440	CoinBigIndex n = `0`;
3441	if (packType) {
3442	double slack = (static_cast<double> (maxSize_-size_-numberAdded))/
3443	static_cast<double> (majorDim_);
3444	slack = CoinMax(`0.0`,slack-`0.01`);
3445	if (!slack) {
3446	for (int i = `0`; i < majorDim_; ++i) {
3447	int thisCount = newStart[i];
3448	newStart[i]=n;
3449	n += length_[i] + thisCount;
3450	}
3451	} else {
3452	double added=`0.0`;
3453	for (int i = `0`; i < majorDim_; ++i) {
3454	int thisCount = newStart[i];
3455	newStart[i]=n;
3456	added += slack;
3457	double extra=`0`;
3458	if (added>=`1.0`) {
3459	extra = floor(added);
3460	added -= extra;
3461	}
3462	n += length_[i] + thisCount+ static_cast<int> (extra);
3463	}
3464	}
3465	newStart[majorDim_]=n;
3466	}
3467	if (packType) {
3468	maxSize_ = CoinMax(maxSize_, n);
3469	int * newIndex = new int[maxSize_];
3470	double * newElem = new double[maxSize_];
3471	for (int i = majorDim_ - `1`; i >= `0`; --i) {
3472	CoinBigIndex start = start_[i];
3473	#ifdef USE_MEMCPY
3474	int length = length_[i];
3475	CoinBigIndex put = newStart[i];
3476	CoinMemcpyN(index_+start,length,newIndex+put);
3477	CoinMemcpyN(element_+start,length,newElem+put);
3478	#else
3479	CoinBigIndex end = start+length_[i];
3480	CoinBigIndex put = newStart[i];
3481	for (CoinBigIndex j=start;j<end;j++) {
3482	newIndex[put]=index_[j];
3483	newElem[put++]=element_[j];
3484	}
3485	#endif
3486	}
3487
3488	delete [] start_;
3489	delete [] index_;
3490	delete [] element_;
3491	start_ = newStart;
3492	index_ = newIndex;
3493	element_ = newElem;
3494	} else if (packType<`0`) {
3495	assert (maxSize_ >= n);
3496	for (int i = majorDim_ - `1`; i >= `0`; --i) {
3497	CoinBigIndex start = start_[i];
3498	int length = length_[i];
3499	CoinBigIndex end = start+length;
3500	CoinBigIndex put = newStart[i];
3501	//if (put==start)
3502	//break;
3503	put += length;
3504	for (CoinBigIndex j=end-`1`;j>=start;j--) {
3505	index_[--put]=index_[j];
3506	element_[put]=element_[j];
3507	}
3508	}
3509	delete [] start_;
3510	start_ = newStart;
3511	} else {
3512	delete[] newStart;
3513	}
3514
3515	// now insert the entries of matrix
3516	for (int i = `0`; i < number; i++) {
3517	CoinBigIndex j;
3518	for ( j=starts[i];j<starts[i+`1`];j++) {
3519	int iIndex = index[j];
3520	element_[start_[iIndex] + length_[iIndex]] = element[j];
3521	index_[start_[iIndex] + (length_[iIndex]++)] = minorDim_;
3522	}
3523	++minorDim_;
3524	}
3525	size_ += starts[number];
3526	#ifndef NDEBUG
3527	int checkSize=`0`;
3528	for (int i=`0`;i<majorDim_;i++) {
3529	checkSize += length_[i];
3530	}
3531	assert (checkSize==size_);
3532	#endif
3533	}
3534
3535	/*
3536	Utility to scan a packed matrix for corruption and inconsistencies. Not
3537	exhaustive, but useful. By default, the method counts coefficients of zero
3538	and reports them, but does not consider them an error. Set zeroesAreError to
3539	true if you want an error.
3540	*/
3541
3542	int CoinPackedMatrix::verifyMtx (int verbosity, bool zeroesAreError) const
3543
3544	{
3545	const double smallCoeff = `1.0e-50` ;
3546	const double largeCoeff = `1.0e50` ;
3547
3548	int majDim = majorDim_ ;
3549	int minDim = minorDim_ ;
3550
3551	std::string majName, minName ;
3552
3553	int m, n ;
3554	if (colOrdered_) {
3555	n = majDim ;
3556	majName = "col" ;
3557	m = minDim ;
3558	minName = "row" ;
3559	} else {
3560	m = majDim ;
3561	majName = "row" ;
3562	n = minDim ;
3563	minName = "col" ;
3564	}
3565
3566	/*
3567	size_ is the number of coefficients, maxSize_ the size of the bulk store.
3568	start_[majDim] should be one past the last valid coefficient in the bulk
3569	store. The actual relation is (#coeffs + #gaps) = start_[majDim].
3570	*/
3571	bool gaps = (size_ < start_[majDim]) ;
3572	CoinBigIndex maxIndex = CoinMin(maxSize_,start_[majDim])-`1` ;
3573
3574	if (verbosity >= `3`) {
3575	std::cout
3576	<< " Matrix is " << ((colOrdered_)?"column":"row") << "-major, "
3577	<< m << " rows X " << n << " cols; " << size_ << " coeffs."
3578	<< std::endl ;
3579	std::cout
3580	<< " Bulk store " << maxSize_ << " coeffs, last coeff at "
3581	<< start_[majDim]-`1` << ", ex maj " << extraMajor_
3582	<< ", ex gap " << extraGap_ ;
3583	if (gaps) std::cout << "; matrix has gaps" ;
3584	std::cout << "." << std::endl ;
3585	}
3586
3587	const CoinBigIndex *const majStarts = start_ ;
3588	const int *const majLens = length_ ;
3589	const int *const minInds = index_ ;
3590	const double *const coeffs = element_ ;
3591	/*
3592	Set up arrays to track use of bulk store entries.
3593	*/
3594	int errs = `0` ;
3595	int zeroes = `0` ;
3596	int refCnt = new* int[maxSize_] ;
3597	CoinZeroN(refCnt,maxSize_) ;
3598	bool inGap = new* bool[maxSize_] ;
3599	CoinZeroN(inGap,maxSize_) ;
3600
3601	for (int majndx = `0` ; majndx < majDim ; majndx++) {
3602	/*
3603	Check that the range of indices for the major vector falls within the bulk
3604	store. If any of these checks fail, it's pointless (and possibly unsafe)
3605	to do more with this vector.
3606
3607	Subtle point: Normally, majStarts[majDim] = maxIndex+1 (one past the
3608	end of the bulk store), and majStarts[k], k < majDim, should be a valid
3609	index. But ... if the last major vector (k = majDim-1) has length 0,
3610	then majStarts[k] = maxIndex. This will propagate back through multiple
3611	major vectors of length 0. Hence the check for length = 0.
3612	*/
3613	CoinBigIndex majStart = majStarts[majndx] ;
3614	int majLen = majLens[majndx] ;
3615
3616	if (majStart < `0` \|\| (majStart == (maxIndex+`1`) && majLen != `0`) \|\|
3617	majStart > maxIndex+`1`) {
3618	if (verbosity >= `1`) {
3619	std::cout
3620	<< " " << majName << " " << majndx
3621	<< ": start " << majStart << " should be between 0 and "
3622	<< maxIndex << "." << std::endl ;
3623	}
3624	errs++ ;
3625	if (majStart >= maxSize_) {
3626	std::cout
3627	<< " " << "index exceeds bulk store limit " << maxSize_
3628	<< "!" << std::endl ;
3629	}
3630	continue ;
3631	}
3632	if (majLen < `0` \|\| majLen > minDim) {
3633	if (verbosity >= `1`) {
3634	std::cout
3635	<< " " << majName << " " << majndx << ": vector length "
3636	<< majLen << " should be between 0 and " << minDim
3637	<< std::endl ;
3638	}
3639	errs++ ;
3640	continue ;
3641	}
3642	CoinBigIndex majEnd = majStart+majLen ;
3643	if (majEnd < `0` \|\| majEnd > maxIndex+`1`) {
3644	if (verbosity >= `1`) {
3645	std::cout
3646	<< " " << majName << " " << majndx
3647	<< ": end " << majEnd << " should be between 0 and "
3648	<< maxIndex << "." << std::endl ;
3649	}
3650	errs++ ;
3651	if (majEnd >= maxSize_) {
3652	std::cout
3653	<< " " << "index exceeds bulk store limit " << maxSize_
3654	<< "!" << std::endl ;
3655	}
3656	continue ;
3657	}
3658	/*
3659	Check that the major vector length is consistent with the distance between
3660	majStart[majndx] and majStart[majndx+1]. If the matrix is gap-free, they
3661	should be equal. We've already confirmed that majStart+majLen is within the
3662	bulk store, so we can continue even if these checks fail.
3663
3664	Recall that the final entry in the major vector start array is one past the
3665	end of the bulk store. The previous tests will check more carefully if
3666	majndx+1 is not the final entry.
3667	*/
3668	CoinBigIndex majStartp1 = majStarts[majndx+`1`] ;
3669	CoinBigIndex startDist = majStartp1-majStart ;
3670	if (majStartp1 < `0` \|\| majStartp1 > maxIndex+`1`) {
3671	if (verbosity >= `1`) {
3672	std::cout
3673	<< " " << majName << " " << majndx
3674	<< ": start of next " << majName << " " << majStartp1
3675	<< " should be between 0 and " << maxIndex+`1` << "." << std::endl ;
3676	}
3677	errs++ ;
3678	if (majStartp1 >= maxSize_) {
3679	std::cout
3680	<< " " << "index exceeds bulk store limit " << maxSize_
3681	<< "!" << std::endl ;
3682	}
3683	} else if ((startDist < `0`) \|\| ((startDist > minDim) && !gaps)) {
3684	if (verbosity >= `1`) {
3685	std::cout
3686	<< " " << majName << " " << majndx << ": distance between "
3687	<< majName << " starts " << startDist
3688	<< " should be between 0 and " << minDim << "." << std::endl ;
3689	}
3690	errs++ ;
3691	} else if (majLen > startDist) {
3692	if (verbosity >= `1`) {
3693	std::cout
3694	<< " " << majName << " " << majndx << ": vector length "
3695	<< majLen << " should not be greater than distance between "
3696	<< majName << " starts " << startDist << std::endl ;
3697	}
3698	errs++ ;
3699	} else if (majLen != startDist && !gaps) {
3700	if (verbosity >= `1`) {
3701	std::cout
3702	<< " " << majName << " " << majndx
3703	<< ": " << majName << " length " << majLen
3704	<< " should equal distance " << startDist << " between "
3705	<< majName << " starts in gap-free matrix." << std::endl ;
3706	}
3707	errs++ ;
3708	}
3709	/*
3710	Scan the major dimension vector, checking for obviously bogus minor indices
3711	and coefficients. Generate reference counts for each bulk store entry.
3712	*/
3713	for (CoinBigIndex ii = majStart ; ii < majEnd ; ii++) {
3714	refCnt[ii]++ ;
3715	int minndx = minInds[ii] ;
3716	if (minndx < `0` \|\| minndx >= minDim) {
3717	if (verbosity >= `1`) {
3718	std::cout
3719	<< " " << majName << " " << majndx << ": "
3720	<< minName << " index " << ii << " is " << minndx
3721	<< ", should be between 0 and " << minDim-`1` << "." << std::endl ;
3722	}
3723	errs++ ;
3724	}
3725	double aij = coeffs[ii] ;
3726	if (CoinIsnan(aij) \|\| CoinAbs(aij) > largeCoeff) {
3727	if (verbosity >= `1`) {
3728	std::cout
3729	<< " (" << ii << ") a<" << majndx << "," << minndx << "> = "
3730	<< aij << " appears bogus." << std::endl ;
3731	}
3732	errs++ ;
3733	}
3734	if (CoinAbs(aij) < smallCoeff) {
3735	if (verbosity >= `4` \|\| zeroesAreError) {
3736	std::cout
3737	<< " (" << ii << ") a<" << majndx << "," << minndx << "> = "
3738	<< aij << " appears bogus." << std::endl ;
3739	}
3740	zeroes++ ;
3741	}
3742	}
3743	/*
3744	And mark the gaps, if any.
3745	*/
3746	if (gaps) {
3747	for (CoinBigIndex ii = majEnd ; ii < majStartp1 ; ii++)
3748	inGap[ii] = true ;
3749	}
3750	}
3751	/*
3752	Check the reference counts. They should all be 1 unless the entry is in a
3753	gap, in which case it should be zero. Anything else is a problem. Allow that
3754	the matrix may not use the full size of the bulk store.
3755	*/
3756	for (CoinBigIndex ii = `0` ; ii <= maxIndex ; ii++) {
3757	if (!((refCnt[ii] == `1` && inGap[ii] == false) \|\|
3758	(refCnt[ii] == `0` && inGap[ii] == true))) {
3759	if (verbosity >= `1`) {
3760	std::cout
3761	<< " Bulk store entry " << ii << " has reference count "
3762	<< refCnt[ii] << "; should be " << ((inGap[ii])?`0`:`1`) << "."
3763	<< std::endl ;
3764	}
3765	errs++ ;
3766	}
3767	}
3768	delete[] refCnt ;
3769	/*
3770	Report the result.
3771	*/
3772	if (zeroesAreError) errs += zeroes ;
3773	if (errs > `0`) {
3774	if (verbosity >= `1`) {
3775	std::cout << " Detected " << errs << " errors in matrix" ;
3776	if (zeroes) std::cout << " (includes " << zeroes << " zeroes)" ;
3777	std::cout << "." << std::endl ;
3778	}
3779	} else {
3780	if (verbosity >= `2`) {
3781	std::cout << " Matrix verified" ;
3782	if (zeroes) std::cout << " (" << zeroes << " zeroes)" ;
3783	std::cout << "." << std::endl ;
3784	}
3785	}
3786
3787	return (errs) ;
3788	}
3789

Browse the source code of OGDF/src/coin/CoinUtils/CoinPackedMatrix.cpp