/* \$Id\$ */ // Copyright (C) 2000, International Business Machines // Corporation and others. All Rights Reserved. // This code is licensed under the terms of the Eclipse Public License (EPL). #include "CoinUtilsConfig.h" #include #include #include #include #include #include #include "CoinPragma.hpp" #include "CoinSort.hpp" #include "CoinHelperFunctions.hpp" #ifndef CLP_NO_VECTOR #include "CoinPackedVectorBase.hpp" #endif #include "CoinFloatEqual.hpp" #include "CoinPackedMatrix.hpp" #if !defined(COIN_COINUTILS_CHECKLEVEL) #define COIN_COINUTILS_CHECKLEVEL 0 #endif //############################################################################# // T must be an integral type (int, CoinBigIndex, etc.) template static inline T CoinLengthWithExtra(T len, double extraGap) { return static_cast(ceil(len * (1 + extraGap))); } //############################################################################# static inline void CoinTestSortedIndexSet(const int num, const int * sorted, const int maxEntry, const char * testingMethod) { if (sorted[0] < 0 || sorted[num-1] >= maxEntry) throw CoinError("bad index", testingMethod, "CoinPackedMatrix"); if (std::adjacent_find(sorted, sorted + num) != sorted + num) throw CoinError("duplicate index", testingMethod, "CoinPackedMatrix"); } //----------------------------------------------------------------------------- static inline int * CoinTestIndexSet(const int numDel, const int * indDel, const int maxEntry, const char * testingMethod) { if (! CoinIsSorted(indDel, indDel + numDel)) { // if not sorted then sort it, test for consistency and return a pointer // to the sorted array int * sorted = new int[numDel]; CoinMemcpyN(indDel, numDel, sorted); std::sort(sorted, sorted + numDel); CoinTestSortedIndexSet(numDel, sorted, maxEntry, testingMethod); return sorted; } // Otherwise it's already sorted, so just test for consistency and return a // 0 pointer. CoinTestSortedIndexSet(numDel, indDel, maxEntry, testingMethod); return 0; } //############################################################################# void CoinPackedMatrix::reserve(const int newMaxMajorDim, const CoinBigIndex newMaxSize, bool create) { if (newMaxMajorDim > maxMajorDim_) { maxMajorDim_ = newMaxMajorDim; int * oldlength = length_; CoinBigIndex * oldstart = start_; length_ = new int[newMaxMajorDim]; start_ = new CoinBigIndex[newMaxMajorDim+1]; start_[0]=0; if (majorDim_ > 0) { CoinMemcpyN(oldlength, majorDim_, length_); CoinMemcpyN(oldstart, majorDim_ + 1, start_); } if (create) { // create empty vectors CoinFillN(length_+majorDim_,maxMajorDim_-majorDim_,0); CoinFillN(start_+majorDim_+1,maxMajorDim_-majorDim_,0); majorDim_=maxMajorDim_; } delete[] oldlength; delete[] oldstart; } if (newMaxSize > maxSize_) { maxSize_ = newMaxSize; int * oldind = index_; double * oldelem = element_; index_ = new int[newMaxSize]; element_ = new double[newMaxSize]; for (int i = majorDim_ - 1; i >= 0; --i) { CoinMemcpyN(oldind+start_[i], length_[i], index_+start_[i]); CoinMemcpyN(oldelem+start_[i], length_[i], element_+start_[i]); } delete[] oldind; delete[] oldelem; } } //----------------------------------------------------------------------------- void CoinPackedMatrix::clear() { majorDim_ = 0; minorDim_ = 0; size_ = 0; } //############################################################################# //############################################################################# void CoinPackedMatrix::setDimensions(int newnumrows, int newnumcols) { const int numrows = getNumRows(); if (newnumrows < 0) newnumrows = numrows; if (newnumrows < numrows) throw CoinError("Bad new rownum (less than current)", "setDimensions", "CoinPackedMatrix"); const int numcols = getNumCols(); if (newnumcols < 0) newnumcols = numcols; if (newnumcols < numcols) throw CoinError("Bad new colnum (less than current)", "setDimensions", "CoinPackedMatrix"); int numplus = 0; if (isColOrdered()) { minorDim_ = newnumrows; numplus = newnumcols - numcols; } else { minorDim_ = newnumcols; numplus = newnumrows - numrows; } if (numplus > 0) { int* lengths = new int[numplus]; CoinZeroN(lengths, numplus); resizeForAddingMajorVectors(numplus, lengths); delete[] lengths; majorDim_ += numplus; //forgot to change majorDim_ } } //----------------------------------------------------------------------------- void CoinPackedMatrix::setExtraGap(const double newGap) { if (newGap < 0) throw CoinError("negative new extra gap", "setExtraGap", "CoinPackedMatrix"); extraGap_ = newGap; } //----------------------------------------------------------------------------- void CoinPackedMatrix::setExtraMajor(const double newMajor) { if (newMajor < 0) throw CoinError("negative new extra major", "setExtraMajor", "CoinPackedMatrix"); extraMajor_ = newMajor; } //############################################################################# #ifndef CLP_NO_VECTOR void CoinPackedMatrix::appendCol(const CoinPackedVectorBase& vec) { if (colOrdered_) appendMajorVector(vec); else appendMinorVector(vec); } #endif //----------------------------------------------------------------------------- void CoinPackedMatrix::appendCol(const int vecsize, const int *vecind, const double *vecelem) { if (colOrdered_) appendMajorVector(vecsize, vecind, vecelem); else appendMinorVector(vecsize, vecind, vecelem); } //----------------------------------------------------------------------------- #ifndef CLP_NO_VECTOR void CoinPackedMatrix::appendCols(const int numcols, const CoinPackedVectorBase * const * cols) { if (colOrdered_) appendMajorVectors(numcols, cols); else appendMinorVectors(numcols, cols); } #endif //----------------------------------------------------------------------------- int CoinPackedMatrix::appendCols(const int numcols, const CoinBigIndex * columnStarts, const int * row, const double * element, int numberRows) { int numberErrors; if (colOrdered_) { numberErrors=appendMajor(numcols, columnStarts, row, element, numberRows); } else { numberErrors=appendMinor(numcols, columnStarts, row, element, numberRows); } return numberErrors; } //----------------------------------------------------------------------------- #ifndef CLP_NO_VECTOR void CoinPackedMatrix::appendRow(const CoinPackedVectorBase& vec) { if (colOrdered_) appendMinorVector(vec); else appendMajorVector(vec); } #endif //----------------------------------------------------------------------------- void CoinPackedMatrix::appendRow(const int vecsize, const int *vecind, const double *vecelem) { if (colOrdered_) appendMinorVector(vecsize, vecind, vecelem); else appendMajorVector(vecsize, vecind, vecelem); } //----------------------------------------------------------------------------- #ifndef CLP_NO_VECTOR void CoinPackedMatrix::appendRows(const int numrows, const CoinPackedVectorBase * const * rows) { if (colOrdered_) { // make sure enough columns if (numrows == 0) return; int i; int maxDim=-1; for (i = numrows - 1; i >= 0; --i) { const int vecsize = rows[i]->getNumElements(); const int* vecind = rows[i]->getIndices(); for (int j = vecsize - 1; j >= 0; --j) maxDim = CoinMax(maxDim,vecind[j]); } maxDim++; if (maxDim>majorDim_) { setDimensions(minorDim_,maxDim); //int nAdd=maxDim-majorDim_; //int * length = new int[nAdd]; //memset(length,0,nAdd*sizeof(int)); //resizeForAddingMajorVectors(nAdd,length); //delete [] length; } appendMinorVectors(numrows, rows); } else { appendMajorVectors(numrows, rows); } } #endif //----------------------------------------------------------------------------- int CoinPackedMatrix::appendRows(const int numrows, const CoinBigIndex * rowStarts, const int * column, const double * element, int numberColumns) { int numberErrors; if (colOrdered_) { numberErrors=appendMinor(numrows, rowStarts, column, element, numberColumns); } else { numberErrors=appendMajor(numrows, rowStarts, column, element, numberColumns); } return numberErrors; } //############################################################################# void CoinPackedMatrix::rightAppendPackedMatrix(const CoinPackedMatrix& matrix) { if (colOrdered_) { if (matrix.colOrdered_) { majorAppendSameOrdered(matrix); } else { majorAppendOrthoOrdered(matrix); } } else { if (matrix.colOrdered_) { minorAppendOrthoOrdered(matrix); } else { minorAppendSameOrdered(matrix); } } } //----------------------------------------------------------------------------- void CoinPackedMatrix::bottomAppendPackedMatrix(const CoinPackedMatrix& matrix) { if (colOrdered_) { if (matrix.colOrdered_) { minorAppendSameOrdered(matrix); } else { minorAppendOrthoOrdered(matrix); } } else { if (matrix.colOrdered_) { majorAppendOrthoOrdered(matrix); } else { majorAppendSameOrdered(matrix); } } } //############################################################################# void CoinPackedMatrix::deleteCols(const int numDel, const int * indDel) { if (numDel) { if (colOrdered_) deleteMajorVectors(numDel, indDel); else deleteMinorVectors(numDel, indDel); } } //----------------------------------------------------------------------------- void CoinPackedMatrix::deleteRows(const int numDel, const int * indDel) { if (numDel) { if (colOrdered_) deleteMinorVectors(numDel, indDel); else deleteMajorVectors(numDel, indDel); } } //############################################################################# /* Replace the elements of a vector. The indices remain the same. At most the number specified will be replaced. The index is between 0 and major dimension of matrix */ void CoinPackedMatrix::replaceVector(const int index, const int numReplace, const double * newElements) { if (index >= 0 && index < majorDim_) { int length = (length_[index] < numReplace) ? length_[index] : numReplace; CoinMemcpyN(newElements, length, element_ + start_[index]); } else { #ifdef COIN_DEBUG throw CoinError("bad index", "replaceVector", "CoinPackedMatrix"); #endif } } /* Modify one element of packed matrix. An element may be added. If the new element is zero it will be deleted unless keepZero true */ void CoinPackedMatrix::modifyCoefficient(int row, int column, double newElement, bool keepZero) { int minorIndex,majorIndex; if (colOrdered_) { majorIndex=column; minorIndex=row; } else { minorIndex=column; majorIndex=row; } if (majorIndex >= 0 && majorIndex < majorDim_) { if (minorIndex >= 0 && minorIndex < minorDim_) { CoinBigIndex j; CoinBigIndex end=start_[majorIndex]+length_[majorIndex];; for (j=start_[majorIndex];j=start_[majorIndex+1]) { int * addedEntries = new int[majorDim_]; memset(addedEntries, 0, majorDim_ * sizeof(int)); addedEntries[majorIndex] = 1; resizeForAddingMinorVectors(addedEntries); delete[] addedEntries; } // So where to insert? We're just going to assume that the entries // in the major vector are in increasing order, so we'll insert the // new entry to the last place we can const CoinBigIndex start = start_[majorIndex]; end = start_[majorIndex]+length_[majorIndex]; // recalculate end for (j = end - 1; j >= start; --j) { if (index_[j] < minorIndex) break; index_[j+1] = index_[j]; element_[j+1] = element_[j]; } ++j; index_[j] = minorIndex; element_[j] = newElement; size_++; length_[majorIndex]++; } } else { #ifdef COIN_DEBUG throw CoinError("bad minor index", "modifyCoefficient", "CoinPackedMatrix"); #endif } } else { #ifdef COIN_DEBUG throw CoinError("bad major index", "modifyCoefficient", "CoinPackedMatrix"); #endif } } /* Return one element of packed matrix. This works for either ordering If it is not present will return 0.0 */ double CoinPackedMatrix::getCoefficient(int row, int column) const { int minorIndex,majorIndex; if (colOrdered_) { majorIndex=column; minorIndex=row; } else { minorIndex=column; majorIndex=row; } double value=0.0; if (majorIndex >= 0 && majorIndex < majorDim_) { if (minorIndex >= 0 && minorIndex < minorDim_) { CoinBigIndex j; CoinBigIndex end=start_[majorIndex]+length_[majorIndex];; for (j=start_[majorIndex];j=threshold) { element_[k]=element_[j]; index_[k++]=index_[j]; } else { eliminatedElement[kbad]=element_[j]; eliminatedIndex[kbad++]=index_[j]; } } if (kbad) { numberEliminated += kbad; length_[i] = k-start_[i]; memcpy(index_+k,eliminatedIndex,kbad*sizeof(int)); memcpy(element_+k,eliminatedElement,kbad*sizeof(double)); } } size_ -= numberEliminated; delete [] eliminatedIndex; delete [] eliminatedElement; return numberEliminated; } //############################################################################# /* Eliminate all elements in matrix whose absolute value is less than threshold.ALSO removes duplicates The column starts are not affected. Returns number of elements eliminated. */ int CoinPackedMatrix::eliminateDuplicates(double threshold) { CoinBigIndex numberEliminated =0; // space for eliminated int * mark = new int [minorDim_]; int i; for (i=0;i=threshold) { element_[k]=element_[j]; index_[k++]=index_[j]; } } numberEliminated += end-k; length_[i] = k-start_[i]; } size_ -= numberEliminated; delete [] mark; return numberEliminated; } //############################################################################# void CoinPackedMatrix::removeGaps(double removeValue) { if (removeValue<0.0) { if (size_size) break; } for (; i < majorDim_; ++i) { const CoinBigIndex si = start_[i]; const int li = length_[i]; start_[i] = size; for (CoinBigIndex j=si;jremoveValue) { index_[put]=index_[j]; element_[put++]=value; } } length_[i]=put-start_[i]; start_[i+1] = put; } size_ = put; } } //############################################################################# /* Really clean up matrix. a) eliminate all duplicate AND small elements in matrix b) remove all gaps and set extraGap_ and extraMajor_ to 0.0 c) reallocate arrays and make max lengths equal to lengths d) orders elements returns number of elements eliminated */ int CoinPackedMatrix::cleanMatrix(double threshold) { if (!majorDim_) { extraGap_=0.0; extraMajor_=0.0; return 0; } CoinBigIndex numberEliminated =0; // space for eliminated int * mark = new int [minorDim_]; int i; for (i=0;i=threshold) { element_[n]=element_[j]; index_[n++]=index_[j]; k++; } } numberEliminated += end-k; length_[i] = n-start_[i]; // sort CoinSort_2(index_+start_[i],index_+n,element_+start_[i]); } start_[majorDim_]=n; size_ -= numberEliminated; assert (n==size_); delete [] mark; extraGap_=0.0; extraMajor_=0.0; maxMajorDim_=majorDim_; maxSize_=size_; // Now reallocate - do smallest ones first int * temp = CoinCopyOfArray(length_,majorDim_); delete [] length_; length_ = temp; CoinBigIndex * temp2 = CoinCopyOfArray(start_,majorDim_+1); delete [] start_; start_ = temp2; temp = CoinCopyOfArray(index_,size_); delete [] index_; index_ = temp; double * temp3 = CoinCopyOfArray(element_,size_); delete [] element_; element_ = temp3; return numberEliminated; } //############################################################################# void CoinPackedMatrix::submatrixOf(const CoinPackedMatrix& matrix, const int numMajor, const int * indMajor) { int i; int* sortedIndPtr = CoinTestIndexSet(numMajor, indMajor, matrix.majorDim_, "submatrixOf"); const int * sortedInd = sortedIndPtr == 0 ? indMajor : sortedIndPtr; gutsOfDestructor(); // Count how many nonzeros there'll be CoinBigIndex nzcnt = 0; const int* length = matrix.getVectorLengths(); for (i = 0; i < numMajor; ++i) { nzcnt += length[sortedInd[i]]; } colOrdered_ = matrix.colOrdered_; maxMajorDim_ = int(numMajor * (1+extraMajor_) + 1); maxSize_ = static_cast (nzcnt * (1+extraMajor_) * (1+extraGap_) + 100); length_ = new int[maxMajorDim_]; start_ = new CoinBigIndex[maxMajorDim_+1]; start_[0]=0; index_ = new int[maxSize_]; element_ = new double[maxSize_]; majorDim_ = 0; minorDim_ = matrix.minorDim_; size_ = 0; #ifdef CLP_NO_VECTOR for (i = 0; i < numMajor; ++i) { int j = sortedInd[i]; CoinBigIndex start = matrix.start_[j]; appendMajorVector(matrix.length_[j],matrix.index_+start,matrix.element_+start); } #else for (i = 0; i < numMajor; ++i) { const CoinShallowPackedVector reqdBySunCC = matrix.getVector(sortedInd[i]) ; appendMajorVector(reqdBySunCC); } #endif delete[] sortedIndPtr; } //############################################################################# void CoinPackedMatrix::submatrixOfWithDuplicates(const CoinPackedMatrix& matrix, const int numMajor, const int * indMajor) { int i; // we allow duplicates - can be useful #ifndef NDEBUG for (i=0; i=matrix.majorDim_) throw CoinError("bad index", "submatrixOfWithDuplicates", "CoinPackedMatrix"); } #endif gutsOfDestructor(); // Get rid of gaps extraMajor_ = 0; extraGap_ = 0; colOrdered_ = matrix.colOrdered_; maxMajorDim_ = numMajor ; const int* length = matrix.getVectorLengths(); length_ = new int[maxMajorDim_]; start_ = new CoinBigIndex[maxMajorDim_+1]; // Count how many nonzeros there'll be CoinBigIndex nzcnt = 0; for (i = 0; i < maxMajorDim_; ++i) { start_[i]=nzcnt; int thisLength = length[indMajor[i]]; nzcnt += thisLength; length_[i]=thisLength; } start_[maxMajorDim_]=nzcnt; maxSize_ = nzcnt ; index_ = new int[maxSize_]; element_ = new double[maxSize_]; majorDim_ = maxMajorDim_; minorDim_ = matrix.minorDim_; size_ = 0; const CoinBigIndex * startOld = matrix.start_; const double * elementOld = matrix.element_; const int * indexOld = matrix.index_; for (i = 0; i < maxMajorDim_; ++i) { int j = indMajor[i]; CoinBigIndex start = startOld[j]; int thisLength = length_[i]; const double * element = elementOld+start; const int * index = indexOld+start; for (int j=0;j=rhs.majorDim_&&maxSize_>=rhs.size_) { majorDim_ = rhs.majorDim_; minorDim_ = rhs.minorDim_; size_ = rhs.size_; extraGap_ = rhs.extraGap_; extraMajor_ = rhs.extraMajor_; CoinMemcpyN(rhs.length_, majorDim_,length_); CoinMemcpyN(rhs.start_, majorDim_+1,start_); if (size_==start_[majorDim_]) { CoinMemcpyN(rhs.index_ , size_, index_); CoinMemcpyN(rhs.element_ , size_, element_); } else { // we can't just simply memcpy these content over, because that can // upset memory debuggers like purify if there were gaps and those gaps // were uninitialized memory blocks for (int i = majorDim_ - 1; i >= 0; --i) { CoinMemcpyN(rhs.index_ + start_[i], length_[i], index_ + start_[i]); CoinMemcpyN(rhs.element_ + start_[i], length_[i], element_ + start_[i]); } } } else { copyOf(rhs); } } //############################################################################# // This method is essentially the same as minorAppendOrthoOrdered(). However, // since we start from an empty matrix, lots of fluff can be avoided. void CoinPackedMatrix::reverseOrderedCopyOf(const CoinPackedMatrix& rhs) { if (this == &rhs) { reverseOrdering(); return; } int i; colOrdered_ = !rhs.colOrdered_; majorDim_ = rhs.minorDim_; minorDim_ = rhs.majorDim_; size_ = rhs.size_; if (size_ == 0) { // we still need to allocate starts and lengths maxMajorDim_=majorDim_; delete[] start_; delete[] length_; delete[] index_; delete[] element_; start_ = new CoinBigIndex[maxMajorDim_ + 1]; length_ = new int[maxMajorDim_]; for (i = 0; i < majorDim_; ++i) { start_[i] = 0; length_[i]=0; } start_[majorDim_]=0; index_ = new int[maxSize_]; element_ = new double[maxSize_]; return; } // Allocate sufficient space (resizeForAddingMinorVectors()) const int newMaxMajorDim_ = CoinMax(maxMajorDim_, CoinLengthWithExtra(majorDim_, extraMajor_)); if (newMaxMajorDim_ > maxMajorDim_) { maxMajorDim_ = newMaxMajorDim_; delete[] start_; delete[] length_; start_ = new CoinBigIndex[maxMajorDim_ + 1]; length_ = new int[maxMajorDim_]; } // first compute how long each major-dimension vector will be int * COIN_RESTRICT orthoLength = length_; rhs.countOrthoLength(orthoLength); start_[0] = 0; if (extraGap_ == 0) { for (i = 0; i < majorDim_; ++i) start_[i+1] = start_[i] + orthoLength[i]; } else { const double eg = extraGap_; for (i = 0; i < majorDim_; ++i) start_[i+1] = start_[i] + CoinLengthWithExtra(orthoLength[i], eg); } const CoinBigIndex newMaxSize = CoinMax(maxSize_, CoinLengthWithExtra(getLastStart(), extraMajor_)); if (newMaxSize > maxSize_) { maxSize_ = newMaxSize; delete[] index_; delete[] element_; index_ = new int[maxSize_]; element_ = new double[maxSize_]; # ifdef ZEROFAULT memset(index_,0,(maxSize_*sizeof(int))) ; memset(element_,0,(maxSize_*sizeof(double))) ; # endif } // now insert the entries of matrix minorDim_ = rhs.majorDim_; const CoinBigIndex * COIN_RESTRICT start = rhs.start_; const int * COIN_RESTRICT index = rhs.index_; const int * COIN_RESTRICT length = rhs.length_; const double * COIN_RESTRICT element = rhs.element_; assert (start[0]==0); CoinBigIndex first = 0; for (i = 0; i < minorDim_; ++i) { CoinBigIndex last = first + length[i]; CoinBigIndex j = first; first = start[i+1]; #if 0 if (((last-j)&1)!=0) { const int ind = index[j]; CoinBigIndex put = start_[ind]; start_[ind] = put +1; element_[put] = element[j]; index_[put] = i; j++; } for (; j != last; j+=2) { const int ind0 = index[j]; CoinBigIndex put0 = start_[ind0]; double value0=element[j]; const int ind1 = index[j+1]; CoinBigIndex put1 = start_[ind1]; double value1=element[j+1]; start_[ind0] = put0 +1; start_[ind1] = put1 +1; element_[put0] = value0; index_[put0] = i; element_[put1] = value1; index_[put1] = i; } #else for (; j != last; ++j) { const int ind = index[j]; CoinBigIndex put = start_[ind]; start_[ind] = put +1; element_[put] = element[j]; index_[put] = i; } #endif } // and re-adjust start_ for (i = 0; i < majorDim_; ++i) { start_[i] -= length_[i]; } } //############################################################################# void CoinPackedMatrix::assignMatrix(const bool colordered, const int minor, const int major, const CoinBigIndex numels, double *& elem, int *& ind, CoinBigIndex *& start, int *& len, const int maxmajor, const CoinBigIndex maxsize) { gutsOfDestructor(); colOrdered_ = colordered; element_ = elem; index_ = ind; start_ = start; majorDim_ = major; minorDim_ = minor; size_ = numels; maxMajorDim_ = maxmajor != -1 ? maxmajor : major; maxSize_ = maxsize != -1 ? maxsize : numels; if (len == NULL) { delete [] length_; length_ = new int[maxMajorDim_]; std::adjacent_difference(start + 1, start + (major + 1), length_); length_[0] -= start[0]; } else { length_ = len; } elem = NULL; ind = NULL; start = NULL; len = NULL; } //############################################################################# CoinPackedMatrix & CoinPackedMatrix::operator=(const CoinPackedMatrix& rhs) { if (this != &rhs) { gutsOfDestructor(); extraGap_=rhs.extraGap_; extraMajor_=rhs.extraMajor_; gutsOfOpEqual(rhs.colOrdered_, rhs.minorDim_, rhs.majorDim_, rhs.size_, rhs.element_, rhs.index_, rhs.start_, rhs.length_); } return *this; } //############################################################################# void CoinPackedMatrix::reverseOrdering() { CoinPackedMatrix m; m.extraGap_ = extraMajor_; m.extraMajor_ = extraGap_; m.reverseOrderedCopyOf(*this); swap(m); } //----------------------------------------------------------------------------- void CoinPackedMatrix::transpose() { colOrdered_ = ! colOrdered_; } //----------------------------------------------------------------------------- void CoinPackedMatrix::swap(CoinPackedMatrix& m) { std::swap(colOrdered_, m.colOrdered_); std::swap(extraGap_, m.extraGap_); std::swap(extraMajor_, m.extraMajor_); std::swap(element_, m.element_); std::swap(index_, m.index_); std::swap(start_, m.start_); std::swap(length_, m.length_); std::swap(majorDim_, m.majorDim_); std::swap(minorDim_, m.minorDim_); std::swap(size_, m.size_); std::swap(maxMajorDim_, m.maxMajorDim_); std::swap(maxSize_, m.maxSize_); } //############################################################################# //############################################################################# void CoinPackedMatrix::times(const double * x, double * y) const { if (colOrdered_) timesMajor(x, y); else timesMinor(x, y); } //----------------------------------------------------------------------------- #ifndef CLP_NO_VECTOR void CoinPackedMatrix::times(const CoinPackedVectorBase& x, double * y) const { if (colOrdered_) timesMajor(x, y); else timesMinor(x, y); } #endif //----------------------------------------------------------------------------- void CoinPackedMatrix::transposeTimes(const double * x, double * y) const { if (colOrdered_) timesMinor(x, y); else timesMajor(x, y); } //----------------------------------------------------------------------------- #ifndef CLP_NO_VECTOR void CoinPackedMatrix::transposeTimes(const CoinPackedVectorBase& x, double * y) const { if (colOrdered_) timesMinor(x, y); else timesMajor(x, y); } #endif //############################################################################# //############################################################################# /* Count the number of entries in every minor-dimension vector and fill in an array containing these lengths. */ void CoinPackedMatrix::countOrthoLength(int * orthoLength) const { CoinZeroN(orthoLength, minorDim_); if (size_!=start_[majorDim_]) { // has gaps for (int i = 0; i =0); ++orthoLength[index_[j]]; } } } else { // no gaps const CoinBigIndex last = start_[majorDim_]; for (CoinBigIndex j = 0; j < last; ++j) { assert( index_[j] < minorDim_ && index_[j]>=0); ++orthoLength[index_[j]]; } } } int * CoinPackedMatrix::countOrthoLength() const { int * orthoLength = new int[minorDim_]; countOrthoLength(orthoLength); return orthoLength; } //############################################################################# /* Returns an array containing major indices. The array is getNumElements long and if getVectorStarts() is 0,2,5 then the array would start 0,0,1,1,1,2... This method is provided to go back from a packed format to a triple format. The returned array is allocated with `new int[]`, free it with `delete[]`. */ int * CoinPackedMatrix::getMajorIndices() const { // Check valid if (!majorDim_||start_[majorDim_]!=size_) return NULL; int * array = new int [size_]; for (int i=0;i(), bind2nd(less(), 0), bind2nd(greater_equal(), minorDim_))) != vecind + vecsize) throw CoinError("out of range index", "appendMajorVector", "CoinPackedMatrix"); #endif #endif if (majorDim_ == maxMajorDim_ || vecsize > maxSize_ - getLastStart()) { resizeForAddingMajorVectors(1, &vecsize); } // got to get this again since it might change! const CoinBigIndex last = getLastStart(); // OK, now just append the major-dimension vector to the end length_[majorDim_] = vecsize; CoinMemcpyN(vecind, vecsize, index_ + last); CoinMemcpyN(vecelem, vecsize, element_ + last); if (majorDim_ == 0) start_[0] = 0; start_[majorDim_ + 1] = CoinMin(last + CoinLengthWithExtra(vecsize, extraGap_), maxSize_ ); // LL: Do we want to allow appending a vector that has more entries than // the current size? if (vecsize > 0) { minorDim_ = CoinMax(minorDim_, (*std::max_element(vecind, vecind+vecsize)) + 1); } ++majorDim_; size_ += vecsize; } //----------------------------------------------------------------------------- #ifndef CLP_NO_VECTOR void CoinPackedMatrix::appendMajorVector(const CoinPackedVectorBase& vec) { appendMajorVector(vec.getNumElements(), vec.getIndices(), vec.getElements()); } //----------------------------------------------------------------------------- void CoinPackedMatrix::appendMajorVectors(const int numvecs, const CoinPackedVectorBase * const * vecs) { int i; CoinBigIndex nz = 0; for (i = 0; i < numvecs; ++i) nz += CoinLengthWithExtra(vecs[i]->getNumElements(), extraGap_); reserve(majorDim_ + numvecs, getLastStart() + nz); for (i = 0; i < numvecs; ++i) appendMajorVector(*vecs[i]); } #endif //############################################################################# void CoinPackedMatrix::appendMinorVector(const int vecsize, const int *vecind, const double *vecelem) { if (vecsize == 0) { ++minorDim_; // empty row/column - still need to increase return; } int i; #if COIN_COINUTILS_CHECKLEVEL > 3 // Test if any of the indices are out of range for (i = 0; i < vecsize; ++i) { if (vecind[i] < 0 || vecind[i] >= majorDim_) throw CoinError("out of range index", "appendMinorVector", "CoinPackedMatrix"); } // Test if there are duplicate indices int* sortedind = CoinCopyOfArray(vecind, vecsize); std::sort(sortedind, sortedind+vecsize); if (std::adjacent_find(sortedind, sortedind+vecsize) != sortedind+vecsize) { throw CoinError("identical indices", "appendMinorVector", "CoinPackedMatrix"); } #endif // test that there's a gap at the end of every major-dimension vector where // we want to add a new entry for (i = vecsize - 1; i >= 0; --i) { const int j = vecind[i]; if (start_[j] + length_[j] == start_[j+1]) break; } if (i >= 0) { int * addedEntries = new int[majorDim_]; memset(addedEntries, 0, majorDim_ * sizeof(int)); for (i = vecsize - 1; i >= 0; --i) addedEntries[vecind[i]] = 1; resizeForAddingMinorVectors(addedEntries); delete[] addedEntries; } // OK, now insert the entries of the minor-dimension vector for (i = vecsize - 1; i >= 0; --i) { const int j = vecind[i]; const CoinBigIndex posj = start_[j] + (length_[j]++); index_[posj] = minorDim_; element_[posj] = vecelem[i]; } ++minorDim_; size_ += vecsize; } //----------------------------------------------------------------------------- #ifndef CLP_NO_VECTOR void CoinPackedMatrix::appendMinorVector(const CoinPackedVectorBase& vec) { appendMinorVector(vec.getNumElements(), vec.getIndices(), vec.getElements()); } //----------------------------------------------------------------------------- void CoinPackedMatrix::appendMinorVectors(const int numvecs, const CoinPackedVectorBase * const * vecs) { if (numvecs == 0) return; int i; int * addedEntries = new int[majorDim_]; CoinZeroN(addedEntries, majorDim_); for (i = numvecs - 1; i >= 0; --i) { const int vecsize = vecs[i]->getNumElements(); const int* vecind = vecs[i]->getIndices(); for (int j = vecsize - 1; j >= 0; --j) { #ifdef COIN_DEBUG if (vecind[j] < 0 || vecind[j] >= majorDim_) throw CoinError("out of range index", "appendMinorVectors", "CoinPackedMatrix"); #endif ++addedEntries[vecind[j]]; } } for (i = majorDim_ - 1; i >= 0; --i) { if (start_[i] + length_[i] + addedEntries[i] > start_[i+1]) break; } if (i >= 0) resizeForAddingMinorVectors(addedEntries); delete[] addedEntries; // now insert the entries of the vectors for (i = 0; i < numvecs; ++i) { const int vecsize = vecs[i]->getNumElements(); const int* vecind = vecs[i]->getIndices(); const double* vecelem = vecs[i]->getElements(); for (int j = vecsize - 1; j >= 0; --j) { const int ind = vecind[j]; element_[start_[ind] + length_[ind]] = vecelem[j]; index_[start_[ind] + (length_[ind]++)] = minorDim_; } ++minorDim_; size_ += vecsize; } } #endif //############################################################################# //############################################################################# void CoinPackedMatrix::majorAppendSameOrdered(const CoinPackedMatrix& matrix) { if (minorDim_ != matrix.minorDim_) { throw CoinError("dimension mismatch", "rightAppendSameOrdered", "CoinPackedMatrix"); } if (matrix.majorDim_ == 0) return; int i; if (majorDim_ + matrix.majorDim_ > maxMajorDim_ || getLastStart() + matrix.getLastStart() > maxSize_) { // we got to resize before we add. note that the resizing method // properly fills out start_ and length_ for the major-dimension // vectors to be added! resizeForAddingMajorVectors(matrix.majorDim_, matrix.length_); start_ += majorDim_; for (i = 0; i < matrix.majorDim_; ++i) { const int l = matrix.length_[i]; CoinMemcpyN(matrix.index_ + matrix.start_[i], l, index_ + start_[i]); CoinMemcpyN(matrix.element_ + matrix.start_[i], l, element_ + start_[i]); } start_ -= majorDim_; } else { start_ += majorDim_; length_ += majorDim_; for (i = 0; i < matrix.majorDim_; ++i) { const int l = matrix.length_[i]; CoinMemcpyN(matrix.index_ + matrix.start_[i], l, index_ + start_[i]); CoinMemcpyN(matrix.element_ + matrix.start_[i], l, element_ + start_[i]); start_[i+1] = start_[i] + matrix.start_[i+1] - matrix.start_[i]; length_[i] = l; } start_ -= majorDim_; length_ -= majorDim_; } majorDim_ += matrix.majorDim_; size_ += matrix.size_; } //----------------------------------------------------------------------------- void CoinPackedMatrix::minorAppendSameOrdered(const CoinPackedMatrix& matrix) { if (majorDim_ != matrix.majorDim_) { throw CoinError("dimension mismatch", "bottomAppendSameOrdered", "CoinPackedMatrix"); } if (matrix.minorDim_ == 0) return; int i; for (i = majorDim_ - 1; i >= 0; --i) { if (start_[i] + length_[i] + matrix.length_[i] > start_[i+1]) break; } if (i >= 0) resizeForAddingMinorVectors(matrix.length_); // now insert the entries of matrix for (i = majorDim_ - 1; i >= 0; --i) { const int l = matrix.length_[i]; std::transform(matrix.index_ + matrix.start_[i], matrix.index_ + (matrix.start_[i] + l), index_ + (start_[i] + length_[i]), std::bind2nd(std::plus(), minorDim_)); CoinMemcpyN(matrix.element_ + matrix.start_[i], l, element_ + (start_[i] + length_[i])); length_[i] += l; } minorDim_ += matrix.minorDim_; size_ += matrix.size_; } //----------------------------------------------------------------------------- void CoinPackedMatrix::majorAppendOrthoOrdered(const CoinPackedMatrix& matrix) { if (minorDim_ != matrix.majorDim_) { throw CoinError("dimension mismatch", "majorAppendOrthoOrdered", "CoinPackedMatrix"); } if (matrix.majorDim_ == 0) return; int i; CoinBigIndex j; // this trickery is needed because MSVC++ is not willing to delete[] a // 'const int *' int * orthoLengthPtr = matrix.countOrthoLength(); const int * orthoLength = orthoLengthPtr; if (majorDim_ + matrix.minorDim_ > maxMajorDim_) { resizeForAddingMajorVectors(matrix.minorDim_, orthoLength); } else { const double extra_gap = extraGap_; start_ += majorDim_; for (i = 0; i < matrix.minorDim_ ; ++i) { start_[i+1] = start_[i] + CoinLengthWithExtra(orthoLength[i], extra_gap); } start_ -= majorDim_; if (start_[majorDim_ + matrix.minorDim_] > maxSize_) { resizeForAddingMajorVectors(matrix.minorDim_, orthoLength); } } // At this point everything is big enough to accommodate the new entries. // Also, start_ is set to the correct starting points for all the new // major-dimension vectors. The length of the new major-dimension vectors // may or may not be correctly set. Hence we just zero them out and they'll // be set when the entries are actually added below. start_ += majorDim_; length_ += majorDim_; CoinZeroN(length_, matrix.minorDim_); for (i = 0; i < matrix.majorDim_; ++i) { const CoinBigIndex last = matrix.getVectorLast(i); for (j = matrix.getVectorFirst(i); j < last; ++j) { const int ind = matrix.index_[j]; element_[start_[ind] + length_[ind]] = matrix.element_[j]; index_[start_[ind] + (length_[ind]++)] = i; } } length_ -= majorDim_; start_ -= majorDim_; // We need to update majorDim_ and size_. We can just add in from matrix majorDim_ += matrix.minorDim_; size_ += matrix.size_; delete[] orthoLengthPtr; } //----------------------------------------------------------------------------- void CoinPackedMatrix::minorAppendOrthoOrdered(const CoinPackedMatrix& matrix) { if (majorDim_ != matrix.minorDim_) { throw CoinError("dimension mismatch", "bottomAppendOrthoOrdered", "CoinPackedMatrix"); } if (matrix.majorDim_ == 0) return; int i; // first compute how many entries will be added to each major-dimension // vector, and if needed, resize the matrix to accommodate all // this trickery is needed because MSVC++ is not willing to delete[] a // 'const int *' int * addedEntriesPtr = matrix.countOrthoLength(); const int * addedEntries = addedEntriesPtr; for (i = majorDim_ - 1; i >= 0; --i) { if (start_[i] + length_[i] + addedEntries[i] > start_[i+1]) break; } if (i >= 0) resizeForAddingMinorVectors(addedEntries); delete[] addedEntriesPtr; // now insert the entries of matrix for (i = 0; i < matrix.majorDim_; ++i) { const CoinBigIndex last = matrix.getVectorLast(i); for (CoinBigIndex j = matrix.getVectorFirst(i); j != last; ++j) { const int ind = matrix.index_[j]; element_[start_[ind] + length_[ind]] = matrix.element_[j]; index_[start_[ind] + (length_[ind]++)] = minorDim_; } ++minorDim_; } size_ += matrix.size_; } //############################################################################# //############################################################################# void CoinPackedMatrix::deleteMajorVectors(const int numDel, const int * indDel) { if (numDel == majorDim_) { // everything is deleted majorDim_ = 0; minorDim_ = 0; size_ = 0; // Get rid of memory as well maxMajorDim_ = 0; delete [] length_; length_ = NULL; delete [] start_; start_ = new CoinBigIndex[1]; start_[0]=0;; delete [] element_; element_=NULL; delete [] index_; index_=NULL; maxSize_ = 0; return; } if (!extraGap_&&!extraMajor_) { // See if this is faster char * keep = new char[majorDim_]; memset(keep,1,majorDim_); for (int i=0;i=0&&k 1) { CoinCopy(start_ + (ind + 1), start_ + ind1, start_ + (ind - i)); CoinCopy(length_ + (ind + 1), length_ + ind1, length_ + (ind - i)); } } // copy the last block of length_ and start_ const int ind = sortedDel[last]; deleted += length_[ind]; if (sortedDel[last] != majorDim_ - 1) { const int ind1 = majorDim_; CoinCopy(start_ + (ind + 1), start_ + ind1, start_ + (ind - last)); CoinCopy(length_ + (ind + 1), length_ + ind1, length_ + (ind - last)); } majorDim_ -= numDel; const int lastlength = CoinLengthWithExtra(length_[majorDim_-1], extraGap_); start_[majorDim_] = CoinMin(start_[majorDim_-1] + lastlength, maxSize_); size_ -= deleted; // if the very first major vector was deleted then copy the new first major // vector to the beginning to make certain that start_[0] is 0. This may // not be necessary, but better safe than sorry... if (sortedDel[0] == 0) { CoinCopyN(index_ + start_[0], length_[0], index_); CoinCopyN(element_ + start_[0], length_[0], element_); start_[0] = 0; } delete[] sortedDelPtr; } } //############################################################################# void CoinPackedMatrix::deleteMinorVectors(const int numDel, const int * indDel) { if (numDel == minorDim_) { // everything is deleted minorDim_ = 0; size_ = 0; // Get rid of as much memory as possible memset(length_,0,majorDim_*sizeof(int)); memset(start_,0,(majorDim_+1)*sizeof(CoinBigIndex )); delete [] element_; element_=NULL; delete [] index_; index_=NULL; maxSize_ = 0; return; } int i, j, k; // first compute the new index of every row int* newindexPtr = new int[minorDim_]; CoinZeroN(newindexPtr, minorDim_); for (j = 0; j < numDel; ++j) { const int ind = indDel[j]; #ifdef COIN_DEBUG if (ind < 0 || ind >= minorDim_) throw CoinError("out of range index", "deleteMinorVectors", "CoinPackedMatrix"); if (newindexPtr[ind] == -1) throw CoinError("duplicate index", "deleteMinorVectors", "CoinPackedMatrix"); #endif newindexPtr[ind] = -1; } for (i = 0, k = 0; i < minorDim_; ++i) { if (newindexPtr[i] != -1) { newindexPtr[i] = k++; } } // Now crawl through the matrix const int * newindex = newindexPtr; #ifdef TAKEOUT int mcount[400]; memset(mcount,0,400*sizeof(int)); for (i = 0; i < majorDim_; ++i) { int * index = index_ + start_[i]; double * elem = element_ + start_[i]; const int length_i = length_[i]; for (j = 0, k = 0; j < length_i; ++j) { mcount[index[j]]++; } } for (i=0;i=0) printf("Keeping original row %d (new %d) with count of %d\n", i,newindex[i],mcount[i]); else printf("deleting row %d with count of %d\n", i,mcount[i]); } } #endif if (!extraGap_) { // pack down size_=0; for (i = 0; i < majorDim_; ++i) { int * index = index_ + start_[i]; double * elem = element_ + start_[i]; start_[i]=size_; const int length_i = length_[i]; for (j = 0; j < length_i; ++j) { const int ind = newindex[index[j]]; if (ind >= 0) { index_[size_] = ind; element_[size_++] = elem[j]; } } length_[i] = size_-start_[i]; } start_[majorDim_]=size_; } else { int deleted = 0; for (i = 0; i < majorDim_; ++i) { int * index = index_ + start_[i]; double * elem = element_ + start_[i]; const int length_i = length_[i]; for (j = 0, k = 0; j < length_i; ++j) { const int ind = newindex[index[j]]; if (ind != -1) { index[k] = ind; elem[k++] = elem[j]; } } deleted += length_i - k; length_[i] = k; } size_ -= deleted; } delete[] newindexPtr; minorDim_ -= numDel; } //############################################################################# //############################################################################# void CoinPackedMatrix::timesMajor(const double * x, double * y) const { memset(y, 0, minorDim_ * sizeof(double)); for (int i = majorDim_ - 1; i >= 0; --i) { const double x_i = x[i]; if (x_i != 0.0) { const CoinBigIndex last = getVectorLast(i); for (CoinBigIndex j = getVectorFirst(i); j < last; ++j) y[index_[j]] += x_i * element_[j]; } } } //----------------------------------------------------------------------------- #ifndef CLP_NO_VECTOR void CoinPackedMatrix::timesMajor(const CoinPackedVectorBase& x, double * y) const { memset(y, 0, minorDim_ * sizeof(double)); for (CoinBigIndex i = x.getNumElements() - 1; i >= 0; --i) { const double x_i = x.getElements()[i]; if (x_i != 0.0) { const int ind = x.getIndices()[i]; const CoinBigIndex last = getVectorLast(ind); for (CoinBigIndex j = getVectorFirst(ind); j < last; ++j) y[index_[j]] += x_i * element_[j]; } } } #endif //----------------------------------------------------------------------------- void CoinPackedMatrix::timesMinor(const double * x, double * y) const { memset(y, 0, majorDim_ * sizeof(double)); for (int i = majorDim_ - 1; i >= 0; --i) { double y_i = 0; const CoinBigIndex last = getVectorLast(i); for (CoinBigIndex j = getVectorFirst(i); j < last; ++j) y_i += x[index_[j]] * element_[j]; y[i] = y_i; } } //----------------------------------------------------------------------------- #ifndef CLP_NO_VECTOR void CoinPackedMatrix::timesMinor(const CoinPackedVectorBase& x, double * y) const { memset(y, 0, majorDim_ * sizeof(double)); for (int i = majorDim_ - 1; i >= 0; --i) { double y_i = 0; const CoinBigIndex last = getVectorLast(i); for (CoinBigIndex j = getVectorFirst(i); j < last; ++j) y_i += x[index_[j]] * element_[j]; y[i] = y_i; } } #endif //############################################################################# //############################################################################# CoinPackedMatrix::CoinPackedMatrix() : colOrdered_(true), extraGap_(0.0), extraMajor_(0.0), element_(0), index_(0), length_(0), majorDim_(0), minorDim_(0), size_(0), maxMajorDim_(0), maxSize_(0) { start_ = new CoinBigIndex[1]; start_[0] = 0; } //----------------------------------------------------------------------------- CoinPackedMatrix::CoinPackedMatrix(const bool colordered, const double extraMajor, const double extraGap) : colOrdered_(colordered), extraGap_(extraGap), extraMajor_(extraMajor), element_(0), index_(0), length_(0), majorDim_(0), minorDim_(0), size_(0), maxMajorDim_(0), maxSize_(0) { start_ = new CoinBigIndex[1]; start_[0] = 0; } //----------------------------------------------------------------------------- CoinPackedMatrix::CoinPackedMatrix(const bool colordered, const int minor, const int major, const CoinBigIndex numels, const double * elem, const int * ind, const CoinBigIndex * start, const int * len, const double extraMajor, const double extraGap) : colOrdered_(colordered), extraGap_(extraGap), extraMajor_(extraMajor), element_(NULL), index_(NULL), start_(NULL), length_(NULL), majorDim_(0), minorDim_(0), size_(0), maxMajorDim_(0), maxSize_(0) { gutsOfOpEqual(colordered, minor, major, numels, elem, ind, start, len); } //----------------------------------------------------------------------------- CoinPackedMatrix::CoinPackedMatrix(const bool colordered, const int minor, const int major, const CoinBigIndex numels, const double * elem, const int * ind, const CoinBigIndex * start, const int * len) : colOrdered_(colordered), extraGap_(0.0), extraMajor_(0.0), element_(NULL), index_(NULL), start_(NULL), length_(NULL), majorDim_(0), minorDim_(0), size_(0), maxMajorDim_(0), maxSize_(0) { gutsOfOpEqual(colordered, minor, major, numels, elem, ind, start, len); } //----------------------------------------------------------------------------- // makes column ordered from triplets and takes out duplicates // will be sorted // // This is an interesting in-place sorting algorithm; // We have triples, and want to sort them so that triples with the same column // are adjacent. // We begin by computing how many entries there are for each column (columnCount) // and using that to compute where each set of column entries will *end* (startColumn). // As we drop entries into place, startColumn is decremented until it contains // the position where the column entries *start*. // The invalid column index -2 means there's a "hole" in that position; // the invalid column index -1 means the entry in that spot is "where it wants to go". // Initially, no one is where they want to go. // Going back to front, // if that entry is where it wants to go // then leave it there // otherwise pick it up (which leaves a hole), and // for as long as you have an entry in your right hand, // - pick up the entry (with your left hand) in the position where the one in // your right hand wants to go; // - pass the entry in your left hand to your right hand; // - was that entry really just the "hole"? If so, stop. // It could be that all the entries get shuffled in the first loop iteration // and all the rest just confirm that everyone is happy where they are. // We never move an entry that is where it wants to go, so entries are moved at // most once. They may not change position if they happen to initially be // where they want to go when the for loop gets to them. // It depends on how many subpermutations the triples initially defined. // Each while loop takes care of one permutation. // The while loop has to stop, because each time around we mark one entry as happy. // We can't run into a happy entry, because we are decrementing the startColumn // all the time, so we must be running into new entries. // Once we've processed all the slots for a column, it cannot be the case that // there are any others that want to go there. // This all means that we eventually must run into the hole. CoinPackedMatrix::CoinPackedMatrix( const bool colordered, const int * indexRow , const int * indexColumn, const double * element, CoinBigIndex numberElements ) : colOrdered_(colordered), extraGap_(0.0), extraMajor_(0.0), element_(NULL), index_(NULL), start_(NULL), length_(NULL), majorDim_(0), minorDim_(0), size_(0), maxMajorDim_(0), maxSize_(0) { CoinAbsFltEq eq; int * colIndices = new int[numberElements]; int * rowIndices = new int[numberElements]; double * elements = new double[numberElements]; CoinCopyN(element,numberElements,elements); if ( colordered ) { CoinCopyN(indexColumn,numberElements,colIndices); CoinCopyN(indexRow,numberElements,rowIndices); } else { CoinCopyN(indexColumn,numberElements,rowIndices); CoinCopyN(indexRow,numberElements,colIndices); } int numberRows; int numberColumns; if (numberElements ) { numberRows = *std::max_element(rowIndices,rowIndices+numberElements)+1; numberColumns = *std::max_element(colIndices,colIndices+numberElements)+1; } else { numberRows = 0; numberColumns = 0; } int * rowCount = new int[numberRows]; int * columnCount = new int[numberColumns]; CoinBigIndex * startColumn = new CoinBigIndex[numberColumns+1]; int * lengths = new int[numberColumns+1]; int iColumn,i; CoinBigIndex k; for (i=0;i=0;k--) { iColumn=colIndices[k]; if (iColumn>=0) { /* pick up the entry with your right hand */ double value = elements[k]; int iRow=rowIndices[k]; colIndices[k]=-2; /* the hole */ while (1) { /* pick this up with your left */ CoinBigIndex iLook=startColumn[iColumn]-1; double valueSave=elements[iLook]; int iColumnSave=colIndices[iLook]; int iRowSave=rowIndices[iLook]; /* put the right-hand entry where it wanted to go */ startColumn[iColumn]=iLook; elements[iLook]=value; rowIndices[iLook]=iRow; colIndices[iLook]=-1; /* mark it as being where it wants to be */ /* there was something there */ if (iColumnSave>=0) { iColumn=iColumnSave; value=valueSave; iRow=iRowSave; } else if (iColumnSave == -2) { /* that was the hole */ break; } else { assert(1==0); /* should never happen */ } /* endif */ } /* endwhile */ } /* endif */ } /* endfor */ /* now pack the elements and combine entries with the same row and column */ /* also, drop entries with "small" coefficients */ numberElements=0; for (iColumn=0;iColumnstart) { int lastRow; double lastValue; // sorts on indices dragging elements with CoinSort_2(rowIndices+start,rowIndices+end,elements+start,CoinFirstLess_2()); lastRow=rowIndices[start]; lastValue=elements[start]; for (i=start+1;ilastRow) { //if(fabs(lastValue)>tolerance) { if(!eq(lastValue,0.0)) { rowIndices[numberElements]=lastRow; elements[numberElements]=lastValue; numberElements++; lengths[iColumn]++; } lastRow=iRow; lastValue=value; } else { lastValue+=value; } /* endif */ } /* endfor */ //if(fabs(lastValue)>tolerance) { if(!eq(lastValue,0.0)) { rowIndices[numberElements]=lastRow; elements[numberElements]=lastValue; numberElements++; lengths[iColumn]++; } } } /* endfor */ startColumn[numberColumns]=numberElements; #if 0 gutsOfOpEqual(colordered,numberRows,numberColumns,numberElements,elements,rowIndices,startColumn,lengths); delete [] rowCount; delete [] columnCount; delete [] startColumn; delete [] lengths; delete [] colIndices; delete [] rowIndices; delete [] elements; #else assignMatrix(colordered,numberRows,numberColumns,numberElements, elements,rowIndices,startColumn,lengths); delete [] rowCount; delete [] columnCount; delete [] lengths; delete [] colIndices; #endif } //----------------------------------------------------------------------------- CoinPackedMatrix::CoinPackedMatrix (const CoinPackedMatrix & rhs) : colOrdered_(true), extraGap_(0.0), extraMajor_(0.0), element_(0), index_(0), start_(0), length_(0), majorDim_(0), minorDim_(0), size_(0), maxMajorDim_(0), maxSize_(0) { bool hasGaps = rhs.size_=0) { maxMajorDim_ = majorDim_+ extraForMajor; maxSize_ = size_ + extraElements; assert (maxMajorDim_>0); assert (maxSize_>0); length_ = new int[maxMajorDim_]; CoinMemcpyN(rhs.length_, majorDim_, length_); start_ = new CoinBigIndex[maxMajorDim_+1]; element_ = new double[maxSize_]; index_ = new int[maxSize_]; bool hasGaps = rhs.size_0) { length_ = new int[maxMajorDim_]; start_ = new CoinBigIndex[maxMajorDim_+1]; if (maxSize_>0) { element_ = new double[maxSize_]; index_ = new int[maxSize_]; } CoinBigIndex size=0; const double * oldElement = rhs.element_; const CoinBigIndex * oldStart = rhs.start_; const int * oldIndex = rhs.index_; const int * oldLength = rhs.length_; CoinBigIndex tooSmallCount=0; for (int i = 0 ; i < majorDim_ ; i++) { start_[i]=size; for (CoinBigIndex j=oldStart[i]; j1.0e-21) { element_[size]=value; index_[size++]=oldIndex[j]; } else { tooSmallCount++; } } length_[i]=size-start_[i]; } start_[majorDim_]=size; assert (size_==size+tooSmallCount); size_ = size; } else { start_ = new CoinBigIndex[1]; start_[0]=0; } } } else { // more complicated colOrdered_ = ! colOrdered_; minorDim_ = rhs.majorDim_; majorDim_ = rhs.minorDim_; maxMajorDim_ = majorDim_ + extraForMajor; maxSize_ = CoinMax(size_ + extraElements,1); assert (maxMajorDim_>0); length_ = new int[maxMajorDim_]; start_ = new CoinBigIndex[maxMajorDim_+1]; element_ = new double[maxSize_]; index_ = new int[maxSize_]; bool hasGaps = rhs.size_=0); ++length_[rhs.index_[j]]; } } } else { // no gaps const CoinBigIndex last = rhs.start_[rhs.majorDim_]; for (CoinBigIndex j = 0; j < last; ++j) { assert( rhs.index_[j] < rhs.minorDim_ && rhs.index_[j]>=0); ++length_[rhs.index_[j]]; } } // Now do starts CoinBigIndex size=0; for (i = 0; i =0 && kRow < rhs.minorDim_) { if (newRow[kRow]<0) { // first time newRow[kRow]=iRow; } else { // duplicate numberDuplicate++; int lastRow = newRow[kRow]; newRow[kRow]=iRow; duplicateRow[iRow] = lastRow; } } else { // bad row numberBad++; } } if (numberBad) throw CoinError("bad minor entries", "subset constructor", "CoinPackedMatrix"); // now get size and check columns size_ = 0; int iColumn; numberBad=0; if (!numberDuplicate) { // No duplicates so can do faster // If not much smaller then use original size if (3*majorDim_>2*rhs.majorDim_&& 3*minorDim_>2*rhs.minorDim_) { // now create arrays maxSize_=CoinMax(static_cast (1),rhs.size_); start_ = new CoinBigIndex [numberColumns+1]; length_ = new int [numberColumns]; index_ = new int[maxSize_]; element_ = new double [maxSize_]; // and fill them size_ = 0; start_[0]=0; for (iColumn=0;iColumn=0 && kColumn =0) { index_[size_] = kRow; element_[size_++] = value; } } } else { // bad column numberBad++; } start_[iColumn+1] = size_; length_[iColumn] = size_ - start_[iColumn]; } if (numberBad) throw CoinError("bad major entries", "subset constructor", "CoinPackedMatrix"); } else { for (iColumn=0;iColumn=0 && kColumn =0) size_++; } } else { // bad column numberBad++; } } if (numberBad) throw CoinError("bad major entries", "subset constructor", "CoinPackedMatrix"); // now create arrays maxSize_=CoinMax(static_cast (1),size_); start_ = new CoinBigIndex [numberColumns+1]; length_ = new int [numberColumns]; index_ = new int[maxSize_]; element_ = new double [maxSize_]; // and fill them size_ = 0; start_[0]=0; for (iColumn=0;iColumn=0) { index_[size_] = kRow; element_[size_++] = value; } } start_[iColumn+1] = size_; length_[iColumn] = size_ - start_[iColumn]; } } } else { for (iColumn=0;iColumn=0 && kColumn =0) { size_++; kRow = duplicateRow[kRow]; } } } else { // bad column numberBad++; } } if (numberBad) throw CoinError("bad major entries", "subset constructor", "CoinPackedMatrix"); // now create arrays maxSize_=CoinMax(static_cast (1),size_); start_ = new CoinBigIndex [numberColumns+1]; length_ = new int [numberColumns]; index_ = new int[maxSize_]; element_ = new double [maxSize_]; // and fill them size_ = 0; start_[0]=0; for (iColumn=0;iColumn=0) { index_[size_] = kRow; element_[size_++] = value; kRow = duplicateRow[kRow]; } } start_[iColumn+1] = size_; length_[iColumn] = size_ - start_[iColumn]; } } delete [] newRow; delete [] duplicateRow; } } //----------------------------------------------------------------------------- CoinPackedMatrix::~CoinPackedMatrix () { gutsOfDestructor(); } //############################################################################# //############################################################################# //############################################################################# void CoinPackedMatrix::gutsOfDestructor() { delete[] length_; delete[] start_; delete[] index_; delete[] element_; length_ = 0; start_ = 0; index_ = 0; element_ = 0; } //############################################################################# void CoinPackedMatrix::gutsOfCopyOf(const bool colordered, const int minor, const int major, const CoinBigIndex numels, const double * elem, const int * ind, const CoinBigIndex * start, const int * len, const double extraMajor, const double extraGap) { colOrdered_ = colordered; majorDim_ = major; minorDim_ = minor; size_ = numels; extraGap_ = extraGap; extraMajor_ = extraMajor; maxMajorDim_ = CoinLengthWithExtra(majorDim_, extraMajor_); if (maxMajorDim_ > 0) { delete [] length_; length_ = new int[maxMajorDim_]; if (len == 0) { std::adjacent_difference(start + 1, start + (major + 1), length_); length_[0] -= start[0]; } else { CoinMemcpyN(len, major, length_); } delete [] start_; start_ = new CoinBigIndex[maxMajorDim_+1]; start_[0]=0; CoinMemcpyN(start, major+1, start_); } else { // empty but be safe delete [] length_; length_ = NULL; delete [] start_; start_ = new CoinBigIndex[1]; start_[0]=0; } maxSize_ = maxMajorDim_ > 0 ? start_[major] : 0; maxSize_ = CoinLengthWithExtra(maxSize_, extraMajor_); if (maxSize_ > 0) { delete [] element_; delete []index_; element_ = new double[maxSize_]; index_ = new int[maxSize_]; // we can't just simply memcpy these content over, because that can // upset memory debuggers like purify if there were gaps and those gaps // were uninitialized memory blocks for (int i = majorDim_ - 1; i >= 0; --i) { CoinMemcpyN(ind + start[i], length_[i], index_ + start_[i]); CoinMemcpyN(elem + start[i], length_[i], element_ + start_[i]); } } } //############################################################################# void CoinPackedMatrix::gutsOfCopyOfNoGaps(const bool colordered, const int minor, const int major, const double * elem, const int * ind, const CoinBigIndex * start) { colOrdered_ = colordered; majorDim_ = major; minorDim_ = minor; size_ = start[majorDim_]; extraGap_ = 0; extraMajor_ = 0; maxMajorDim_ = majorDim_; // delete all arrays delete [] length_; delete [] start_; delete [] element_; delete [] index_; if (maxMajorDim_ > 0) { length_ = new int[maxMajorDim_]; assert (!start[0]); start_ = new CoinBigIndex[maxMajorDim_+1]; start_[0]=0; CoinBigIndex last = 0; for (int i=0;i 0) { element_ = new double[maxSize_]; index_ = new int[maxSize_]; CoinMemcpyN(ind , maxSize_, index_); CoinMemcpyN(elem , maxSize_, element_); } else { element_ = NULL; index_ = NULL; } } //############################################################################# void CoinPackedMatrix::gutsOfOpEqual(const bool colordered, const int minor, const int major, const CoinBigIndex numels, const double * elem, const int * ind, const CoinBigIndex * start, const int * len) { colOrdered_ = colordered; majorDim_ = major; minorDim_ = minor; size_ = numels; if (!len && numels > 0 && numels==start[major] && start[0]==0) { // No gaps - do faster if (major>maxMajorDim_||!start_) { maxMajorDim_ = major; delete [] length_; length_ = new int[maxMajorDim_]; delete [] start_; start_ = new CoinBigIndex[maxMajorDim_+1]; } CoinMemcpyN(start,major+1,start_); std::adjacent_difference(start + 1, start + (major + 1), length_); if (numels>maxSize_||!element_) { maxSize_=numels; delete [] element_; delete [] index_; element_ = new double[maxSize_]; index_ = new int[maxSize_]; } CoinMemcpyN(ind,numels,index_); CoinMemcpyN(elem,numels,element_); } else { maxMajorDim_ = CoinLengthWithExtra(majorDim_, extraMajor_); int i; if (maxMajorDim_ > 0) { delete [] length_; length_ = new int[maxMajorDim_]; if (len == 0) { std::adjacent_difference(start + 1, start + (major + 1), length_); length_[0] -= start[0]; } else { CoinMemcpyN(len, major, length_); } delete [] start_; start_ = new CoinBigIndex[maxMajorDim_+1]; start_[0] = 0; if (extraGap_ == 0) { for (i = 0; i < major; ++i) start_[i+1] = start_[i] + length_[i]; } else { const double extra_gap = extraGap_; for (i = 0; i < major; ++i) start_[i+1] = start_[i] + CoinLengthWithExtra(length_[i], extra_gap); } } else { // empty matrix delete [] start_; start_ = new CoinBigIndex[1]; start_[0] = 0; } maxSize_ = maxMajorDim_ > 0 ? start_[major] : 0; maxSize_ = CoinLengthWithExtra(maxSize_, extraMajor_); if (maxSize_ > 0) { delete [] element_; delete [] index_; element_ = new double[maxSize_]; index_ = new int[maxSize_]; assert (maxSize_>=start_[majorDim_-1]+length_[majorDim_-1]); // we can't just simply memcpy these content over, because that can // upset memory debuggers like purify if there were gaps and those gaps // were uninitialized memory blocks for (i = majorDim_ - 1; i >= 0; --i) { CoinMemcpyN(ind + start[i], length_[i], index_ + start_[i]); CoinMemcpyN(elem + start[i], length_[i], element_ + start_[i]); } } } #ifndef NDEBUG for (int i = majorDim_ - 1; i >= 0; --i) { const CoinBigIndex last = getVectorLast(i); for (CoinBigIndex j = getVectorFirst(i); j < last; ++j) { int index = index_[j]; assert (index>=0&&index= 0; --i) { CoinMemcpyN(index_ + start_[i], length_[i], newIndex + newStart[i]); CoinMemcpyN(element_ + start_[i], length_[i], newElem + newStart[i]); } gutsOfDestructor(); start_ = newStart; length_ = newLength; index_ = newIndex; element_ = newElem; } //############################################################################# void CoinPackedMatrix::resizeForAddingMinorVectors(const int * addedEntries) { int i; maxMajorDim_ = CoinMax(CoinLengthWithExtra(majorDim_, extraMajor_), maxMajorDim_); CoinBigIndex * newStart = new CoinBigIndex[maxMajorDim_ + 1]; int * newLength = new int[maxMajorDim_]; // increase the lengths temporarily so that the correct new start positions // can be easily computed (it's faster to modify the lengths and reset them // than do a test for every entry when the start positions are computed. for (i = majorDim_ - 1; i >= 0; --i) newLength[i] = length_[i] + addedEntries[i]; newStart[0] = 0; if (extraGap_ == 0) { for (i = 0; i < majorDim_; ++i) newStart[i+1] = newStart[i] + newLength[i]; } else { const double eg = extraGap_; for (i = 0; i < majorDim_; ++i) newStart[i+1] = newStart[i] + CoinLengthWithExtra(newLength[i], eg); } // reset the lengths for (i = majorDim_ - 1; i >= 0; --i) newLength[i] -= addedEntries[i]; maxSize_ = CoinMax(maxSize_, CoinLengthWithExtra(newStart[majorDim_], extraMajor_)); int * newIndex = new int[maxSize_]; double * newElem = new double[maxSize_]; for (i = majorDim_ - 1; i >= 0; --i) { CoinMemcpyN(index_ + start_[i], length_[i], newIndex + newStart[i]); CoinMemcpyN(element_ + start_[i], length_[i], newElem + newStart[i]); } gutsOfDestructor(); start_ = newStart; length_ = newLength; index_ = newIndex; element_ = newElem; } //############################################################################# //############################################################################# void CoinPackedMatrix::dumpMatrix(const char* fname) const { if (! fname) { printf("Dumping matrix...\n\n"); printf("colordered: %i\n", isColOrdered() ? 1 : 0); const int major = getMajorDim(); const int minor = getMinorDim(); printf("major: %i minor: %i\n", major, minor); for (int i = 0; i < major; ++i) { printf("vec %i has length %i with entries:\n", i, length_[i]); for (CoinBigIndex j = start_[i]; j < start_[i] + length_[i]; ++j) { printf(" %15i %40.25f\n", index_[j], element_[j]); } } printf("\nFinished dumping matrix\n"); } else { FILE* out = fopen(fname, "w"); fprintf(out, "Dumping matrix...\n\n"); fprintf(out, "colordered: %i\n", isColOrdered() ? 1 : 0); const int major = getMajorDim(); const int minor = getMinorDim(); fprintf(out, "major: %i minor: %i\n", major, minor); for (int i = 0; i < major; ++i) { fprintf(out, "vec %i has length %i with entries:\n", i, length_[i]); for (CoinBigIndex j = start_[i]; j < start_[i] + length_[i]; ++j) { fprintf(out, " %15i %40.25f\n", index_[j], element_[j]); } } fprintf(out, "\nFinished dumping matrix\n"); fclose(out); } } void CoinPackedMatrix::printMatrixElement (const int row_val, const int col_val) const { int major_index, minor_index; if (isColOrdered()) { major_index = col_val; minor_index = row_val; } else { major_index = row_val; minor_index = col_val; } if (major_index < 0 || major_index > getMajorDim()-1) { std::cout << "Major index " << major_index << " not in range 0.." << getMajorDim()-1 << std::endl ; } else if (minor_index < 0 || minor_index > getMinorDim()-1) { std::cout << "Minor index " << minor_index << " not in range 0.." << getMinorDim()-1 << std::endl ; } else { CoinBigIndex curr_point = start_[major_index]; const CoinBigIndex stop_point = curr_point+length_[major_index]; double aij = 0.0 ; for ( ; curr_point < stop_point ; curr_point++) { if (index_[curr_point] == minor_index) { aij = element_[curr_point]; break; } } std::cout << aij ; } } #ifndef CLP_NO_VECTOR bool CoinPackedMatrix::isEquivalent2(const CoinPackedMatrix& rhs) const { CoinRelFltEq eq; // Both must be column order or both row ordered and must be of same size if (isColOrdered() ^ rhs.isColOrdered()) { std::cerr<<"Ordering "<= 0; --i) { CoinShallowPackedVector pv = getVector(i); CoinShallowPackedVector rhsPv = rhs.getVector(i); if ( !pv.isEquivalent(rhsPv,eq) ) { std::cerr<<"vector # "< (elems+j); printf("%x %x",xx[0],xx[1]); xx = reinterpret_cast (elems2+j); printf(" %x %x\n",xx[0],xx[1]); } } //return false; } } return true; } #else /* Equivalence. Two matrices are equivalent if they are both by rows or both by columns, they have the same dimensions, and each vector is equivalent. In this method the FloatEqual function operator can be specified. */ bool CoinPackedMatrix::isEquivalent(const CoinPackedMatrix& rhs, const CoinRelFltEq& eq) const { // Both must be column order or both row ordered and must be of same size if ((isColOrdered() ^ rhs.isColOrdered()) || (getNumCols() != rhs.getNumCols()) || (getNumRows() != rhs.getNumRows()) || (getNumElements() != rhs.getNumElements())) return false; const int major = getMajorDim(); const int minor = getMinorDim(); double * values = new double[minor]; memset(values,0,minor*sizeof(double)); bool same=true; for (int i = 0; i < major; ++i) { int length = length_[i]; if (length!=rhs.length_[i]) { same=false; break; } else { CoinBigIndex j; for ( j = start_[i]; j < start_[i] + length; ++j) { int index = index_[j]; values[index]=element_[j]; } for ( j = rhs.start_[i]; j < rhs.start_[i] + length; ++j) { int index = index_[j]; double oldValue = values[index]; values[index]=0.0; if (!eq(oldValue,rhs.element_[j])) { same=false; break; } } if (!same) break; } } delete [] values; return same; } #endif bool CoinPackedMatrix::isEquivalent(const CoinPackedMatrix& rhs) const { return isEquivalent(rhs,CoinRelFltEq()); } /* Sort all columns so indices are increasing.in each column */ void CoinPackedMatrix::orderMatrix() { for (int i=0;i0) or duplicates This version is easy one i.e. adding columns to column ordered */ int CoinPackedMatrix::appendMajor(const int number, const CoinBigIndex * starts, const int * index, const double * element, int numberOther) { int i; int numberErrors=0; CoinBigIndex numberElements = starts[number]; if (majorDim_ + number > maxMajorDim_ || getLastStart() + numberElements > maxSize_) { // we got to resize before we add. note that the resizing method // properly fills out start_ and length_ for the major-dimension // vectors to be added! if (!extraGap_&&!extraMajor_&&numberOther<=0&&!hasGaps()) { // can do faster if (majorDim_+number>maxMajorDim_) { maxMajorDim_ = majorDim_+number; int * newLength = new int[maxMajorDim_]; CoinMemcpyN(length_, majorDim_, newLength); delete [] length_; length_ = newLength; CoinBigIndex * newStart = new CoinBigIndex[maxMajorDim_ + 1]; CoinMemcpyN(start_, majorDim_+1, newStart); delete [] start_; start_ = newStart; } if (size_+numberElements>maxSize_) { maxSize_ = size_+numberElements; double * newElem = new double[maxSize_]; CoinMemcpyN(element_,size_,newElem); delete [] element_; element_ = newElem; int * newIndex = new int[maxSize_]; CoinMemcpyN(index_,size_,newIndex); delete [] index_; index_ = newIndex; } CoinMemcpyN(index,numberElements,index_+size_); // Do minor dimension int lastMinor=-1; for (CoinBigIndex j=0;j0) { char * which = new char[numberOther]; memset(which,0,numberOther); for (i = 0; i < number; i++) { CoinBigIndex put = start_[majorDim_+i]; CoinBigIndex j; for ( j=starts[i];j=0&&iIndex=0&&iIndex0) { char * which = new char[numberOther]; memset(which,0,numberOther); for (i = 0; i < number; i++) { CoinBigIndex put = start_[majorDim_+i]; CoinBigIndex j; for ( j=starts[i];j=0&&iIndex=0&&iIndex0) or duplicates This version is harder one i.e. adding columns to row ordered */ int CoinPackedMatrix::appendMinor(const int number, const CoinBigIndex * starts, const int * index, const double * element, int numberOther) { int i; int numberErrors=0; // first compute how many entries will be added to each major-dimension // vector, and if needed, resize the matrix to accommodate all int * addedEntries = NULL; if (numberOther>0) { addedEntries = new int[majorDim_]; CoinZeroN(addedEntries,majorDim_); numberOther=majorDim_; char * which = new char[numberOther]; memset(which,0,numberOther); for (i = 0; i < number; i++) { CoinBigIndex j; for ( j=starts[i];j=0&&iIndex=0&&iIndexmajorDim_) { if (isColOrdered()) setDimensions(-1,largest+1); else setDimensions(largest+1,-1); } addedEntries = new int[majorDim_]; CoinZeroN(addedEntries,majorDim_); // no checking for (i = 0; i < number; i++) { CoinBigIndex j; for ( j=starts[i];j= 0; i--) { if (start_[i] + length_[i] + addedEntries[i] > start_[i+1]) break; } if (i >= 0) resizeForAddingMinorVectors(addedEntries); delete[] addedEntries; // now insert the entries of matrix for (i = 0; i < number; i++) { CoinBigIndex j; for ( j=starts[i];j= 0; i--) { CoinBigIndex start = start_[i]; if (start + length_[i] + newStart[i] <= nextStart) { nextStart=start; } else { packType=-1; #ifdef ADD_ROW_ANALYZE nBad++; #else break; #endif } } } else { // Need more space packType=1; } #ifdef ADD_ROW_ANALYZE if (!hasGaps()) xxxxxx[9]++; if (packType==-1&&nBad<6) packType=nBad+1; xxxxxx[packType+2]++; if ((xxxxxx[0]%100)==0) { printf("Append "); for (int i=0;i<10;i++) printf("%d ",xxxxxx[i]); printf("\n"); } #endif if (hasGaps()&&packType) packType=1; CoinBigIndex n = 0; if (packType) { double slack = (static_cast (maxSize_-size_-numberAdded))/ static_cast (majorDim_); slack = CoinMax(0.0,slack-0.01); if (!slack) { for (int i = 0; i < majorDim_; ++i) { int thisCount = newStart[i]; newStart[i]=n; n += length_[i] + thisCount; } } else { double added=0.0; for (int i = 0; i < majorDim_; ++i) { int thisCount = newStart[i]; newStart[i]=n; added += slack; double extra=0; if (added>=1.0) { extra = floor(added); added -= extra; } n += length_[i] + thisCount+ static_cast (extra); } } newStart[majorDim_]=n; } if (packType) { maxSize_ = CoinMax(maxSize_, n); int * newIndex = new int[maxSize_]; double * newElem = new double[maxSize_]; for (int i = majorDim_ - 1; i >= 0; --i) { CoinBigIndex start = start_[i]; #ifdef USE_MEMCPY int length = length_[i]; CoinBigIndex put = newStart[i]; CoinMemcpyN(index_+start,length,newIndex+put); CoinMemcpyN(element_+start,length,newElem+put); #else CoinBigIndex end = start+length_[i]; CoinBigIndex put = newStart[i]; for (CoinBigIndex j=start;j= n); for (int i = majorDim_ - 1; i >= 0; --i) { CoinBigIndex start = start_[i]; int length = length_[i]; CoinBigIndex end = start+length; CoinBigIndex put = newStart[i]; //if (put==start) //break; put += length; for (CoinBigIndex j=end-1;j>=start;j--) { index_[--put]=index_[j]; element_[put]=element_[j]; } } delete [] start_; start_ = newStart; } else { delete[] newStart; } // now insert the entries of matrix for (int i = 0; i < number; i++) { CoinBigIndex j; for ( j=starts[i];j= 3) { std::cout << " Matrix is " << ((colOrdered_)?"column":"row") << "-major, " << m << " rows X " << n << " cols; " << size_ << " coeffs." << std::endl ; std::cout << " Bulk store " << maxSize_ << " coeffs, last coeff at " << start_[majDim]-1 << ", ex maj " << extraMajor_ << ", ex gap " << extraGap_ ; if (gaps) std::cout << "; matrix has gaps" ; std::cout << "." << std::endl ; } const CoinBigIndex *const majStarts = start_ ; const int *const majLens = length_ ; const int *const minInds = index_ ; const double *const coeffs = element_ ; /* Set up arrays to track use of bulk store entries. */ int errs = 0 ; int zeroes = 0 ; int *refCnt = new int[maxSize_] ; CoinZeroN(refCnt,maxSize_) ; bool *inGap = new bool[maxSize_] ; CoinZeroN(inGap,maxSize_) ; for (int majndx = 0 ; majndx < majDim ; majndx++) { /* Check that the range of indices for the major vector falls within the bulk store. If any of these checks fail, it's pointless (and possibly unsafe) to do more with this vector. Subtle point: Normally, majStarts[majDim] = maxIndex+1 (one past the end of the bulk store), and majStarts[k], k < majDim, should be a valid index. But ... if the last major vector (k = majDim-1) has length 0, then majStarts[k] = maxIndex. This will propagate back through multiple major vectors of length 0. Hence the check for length = 0. */ CoinBigIndex majStart = majStarts[majndx] ; int majLen = majLens[majndx] ; if (majStart < 0 || (majStart == (maxIndex+1) && majLen != 0) || majStart > maxIndex+1) { if (verbosity >= 1) { std::cout << " " << majName << " " << majndx << ": start " << majStart << " should be between 0 and " << maxIndex << "." << std::endl ; } errs++ ; if (majStart >= maxSize_) { std::cout << " " << "index exceeds bulk store limit " << maxSize_ << "!" << std::endl ; } continue ; } if (majLen < 0 || majLen > minDim) { if (verbosity >= 1) { std::cout << " " << majName << " " << majndx << ": vector length " << majLen << " should be between 0 and " << minDim << std::endl ; } errs++ ; continue ; } CoinBigIndex majEnd = majStart+majLen ; if (majEnd < 0 || majEnd > maxIndex+1) { if (verbosity >= 1) { std::cout << " " << majName << " " << majndx << ": end " << majEnd << " should be between 0 and " << maxIndex << "." << std::endl ; } errs++ ; if (majEnd >= maxSize_) { std::cout << " " << "index exceeds bulk store limit " << maxSize_ << "!" << std::endl ; } continue ; } /* Check that the major vector length is consistent with the distance between majStart[majndx] and majStart[majndx+1]. If the matrix is gap-free, they should be equal. We've already confirmed that majStart+majLen is within the bulk store, so we can continue even if these checks fail. Recall that the final entry in the major vector start array is one past the end of the bulk store. The previous tests will check more carefully if majndx+1 is not the final entry. */ CoinBigIndex majStartp1 = majStarts[majndx+1] ; CoinBigIndex startDist = majStartp1-majStart ; if (majStartp1 < 0 || majStartp1 > maxIndex+1) { if (verbosity >= 1) { std::cout << " " << majName << " " << majndx << ": start of next " << majName << " " << majStartp1 << " should be between 0 and " << maxIndex+1 << "." << std::endl ; } errs++ ; if (majStartp1 >= maxSize_) { std::cout << " " << "index exceeds bulk store limit " << maxSize_ << "!" << std::endl ; } } else if ((startDist < 0) || ((startDist > minDim) && !gaps)) { if (verbosity >= 1) { std::cout << " " << majName << " " << majndx << ": distance between " << majName << " starts " << startDist << " should be between 0 and " << minDim << "." << std::endl ; } errs++ ; } else if (majLen > startDist) { if (verbosity >= 1) { std::cout << " " << majName << " " << majndx << ": vector length " << majLen << " should not be greater than distance between " << majName << " starts " << startDist << std::endl ; } errs++ ; } else if (majLen != startDist && !gaps) { if (verbosity >= 1) { std::cout << " " << majName << " " << majndx << ": " << majName << " length " << majLen << " should equal distance " << startDist << " between " << majName << " starts in gap-free matrix." << std::endl ; } errs++ ; } /* Scan the major dimension vector, checking for obviously bogus minor indices and coefficients. Generate reference counts for each bulk store entry. */ for (CoinBigIndex ii = majStart ; ii < majEnd ; ii++) { refCnt[ii]++ ; int minndx = minInds[ii] ; if (minndx < 0 || minndx >= minDim) { if (verbosity >= 1) { std::cout << " " << majName << " " << majndx << ": " << minName << " index " << ii << " is " << minndx << ", should be between 0 and " << minDim-1 << "." << std::endl ; } errs++ ; } double aij = coeffs[ii] ; if (CoinIsnan(aij) || CoinAbs(aij) > largeCoeff) { if (verbosity >= 1) { std::cout << " (" << ii << ") a<" << majndx << "," << minndx << "> = " << aij << " appears bogus." << std::endl ; } errs++ ; } if (CoinAbs(aij) < smallCoeff) { if (verbosity >= 4 || zeroesAreError) { std::cout << " (" << ii << ") a<" << majndx << "," << minndx << "> = " << aij << " appears bogus." << std::endl ; } zeroes++ ; } } /* And mark the gaps, if any. */ if (gaps) { for (CoinBigIndex ii = majEnd ; ii < majStartp1 ; ii++) inGap[ii] = true ; } } /* Check the reference counts. They should all be 1 unless the entry is in a gap, in which case it should be zero. Anything else is a problem. Allow that the matrix may not use the full size of the bulk store. */ for (CoinBigIndex ii = 0 ; ii <= maxIndex ; ii++) { if (!((refCnt[ii] == 1 && inGap[ii] == false) || (refCnt[ii] == 0 && inGap[ii] == true))) { if (verbosity >= 1) { std::cout << " Bulk store entry " << ii << " has reference count " << refCnt[ii] << "; should be " << ((inGap[ii])?0:1) << "." << std::endl ; } errs++ ; } } delete[] refCnt ; /* Report the result. */ if (zeroesAreError) errs += zeroes ; if (errs > 0) { if (verbosity >= 1) { std::cout << " Detected " << errs << " errors in matrix" ; if (zeroes) std::cout << " (includes " << zeroes << " zeroes)" ; std::cout << "." << std::endl ; } } else { if (verbosity >= 2) { std::cout << " Matrix verified" ; if (zeroes) std::cout << " (" << zeroes << " zeroes)" ; std::cout << "." << std::endl ; } } return (errs) ; }