00001
00002
00003
00004 #if defined(_MSC_VER)
00005
00006 # pragma warning(disable:4786)
00007 #endif
00008
00009 #include <algorithm>
00010 #include <numeric>
00011 #include <cassert>
00012 #include <cstdio>
00013 #include <iostream>
00014
00015 #include "CoinSort.hpp"
00016 #include "CoinHelperFunctions.hpp"
00017 #include "CoinPackedVectorBase.hpp"
00018 #include "CoinPackedMatrix.hpp"
00019 #include <stdio.h>
00020
00021
00022
00023 static inline int
00024 CoinLengthWithExtra(int len, double extraGap)
00025 {
00026 return static_cast<int>(ceil(len * (1 + extraGap)));
00027 }
00028
00029
00030
00031 static inline void
00032 CoinTestSortedIndexSet(const int num, const int * sorted, const int maxEntry,
00033 const char * testingMethod)
00034 {
00035 if (sorted[0] < 0 || sorted[num-1] >= maxEntry)
00036 throw CoinError("bad index", testingMethod, "CoinPackedMatrix");
00037 if (std::adjacent_find(sorted, sorted + num) != sorted + num)
00038 throw CoinError("duplicate index", testingMethod, "CoinPackedMatrix");
00039 }
00040
00041
00042
00043 static inline int *
00044 CoinTestIndexSet(const int numDel, const int * indDel, const int maxEntry,
00045 const char * testingMethod)
00046 {
00047 if (! CoinIsSorted(indDel, indDel + numDel)) {
00048
00049
00050 int * sorted = new int[numDel];
00051 CoinMemcpyN(indDel, numDel, sorted);
00052 std::sort(sorted, sorted + numDel);
00053 CoinTestSortedIndexSet(numDel, sorted, maxEntry, testingMethod);
00054 return sorted;
00055 }
00056
00057
00058
00059 CoinTestSortedIndexSet(numDel, indDel, maxEntry, testingMethod);
00060 return 0;
00061 }
00062
00063
00064
00065 void
00066 CoinPackedMatrix::reserve(const int newMaxMajorDim, const CoinBigIndex newMaxSize)
00067 {
00068 if (newMaxMajorDim > maxMajorDim_) {
00069 maxMajorDim_ = newMaxMajorDim;
00070 int * oldlength = length_;
00071 CoinBigIndex * oldstart = start_;
00072 length_ = new int[newMaxMajorDim];
00073 start_ = new CoinBigIndex[newMaxMajorDim+1];
00074 if (majorDim_ > 0) {
00075 CoinMemcpyN(oldlength, majorDim_, length_);
00076 CoinMemcpyN(oldstart, majorDim_ + 1, start_);
00077 }
00078 delete[] oldlength;
00079 delete[] oldstart;
00080 }
00081 if (newMaxSize > maxSize_) {
00082 maxSize_ = newMaxSize;
00083 int * oldind = index_;
00084 double * oldelem = element_;
00085 index_ = new int[newMaxSize];
00086 element_ = new double[newMaxSize];
00087 for (int i = majorDim_ - 1; i >= 0; --i) {
00088 CoinMemcpyN(oldind+start_[i], length_[i], index_+start_[i]);
00089 CoinMemcpyN(oldelem+start_[i], length_[i], element_+start_[i]);
00090 }
00091 delete[] oldind;
00092 delete[] oldelem;
00093 }
00094 }
00095
00096
00097
00098 void
00099 CoinPackedMatrix::clear()
00100 {
00101 majorDim_ = 0;
00102 minorDim_ = 0;
00103 size_ = 0;
00104 }
00105
00106
00107
00108
00109 void
00110 CoinPackedMatrix::setDimensions(int newnumrows, int newnumcols)
00111 throw(CoinError)
00112 {
00113 const int numrows = getNumRows();
00114 if (newnumrows < 0)
00115 newnumrows = numrows;
00116 if (newnumrows < numrows)
00117 throw CoinError("Bad new rownum (less than current)",
00118 "setDimensions", "CoinPackedMatrix");
00119
00120 const int numcols = getNumCols();
00121 if (newnumcols < 0)
00122 newnumcols = numcols;
00123 if (newnumcols < numcols)
00124 throw CoinError("Bad new colnum (less than current)",
00125 "setDimensions", "CoinPackedMatrix");
00126
00127 int numplus = 0;
00128 if (isColOrdered()) {
00129 minorDim_ = newnumrows;
00130 numplus = newnumcols - numcols;
00131 } else {
00132 minorDim_ = newnumcols;
00133 numplus = newnumrows - numrows;
00134 }
00135 if (numplus > 0) {
00136 int* lengths = new int[numplus];
00137 CoinZeroN(lengths, numplus);
00138 resizeForAddingMajorVectors(numplus, lengths);
00139 delete[] lengths;
00140 majorDim_ += numplus;
00141 }
00142
00143 }
00144
00145
00146
00147 void
00148 CoinPackedMatrix::setExtraGap(const double newGap) throw(CoinError)
00149 {
00150 if (newGap < 0)
00151 throw CoinError("negative new extra gap",
00152 "setExtraGap", "CoinPackedMatrix");
00153 extraGap_ = newGap;
00154 }
00155
00156
00157
00158 void
00159 CoinPackedMatrix::setExtraMajor(const double newMajor) throw(CoinError)
00160 {
00161 if (newMajor < 0)
00162 throw CoinError("negative new extra major",
00163 "setExtraMajor", "CoinPackedMatrix");
00164 extraMajor_ = newMajor;
00165 }
00166
00167
00168
00169 void
00170 CoinPackedMatrix::appendCol(const CoinPackedVectorBase& vec) throw(CoinError)
00171 {
00172 if (colOrdered_)
00173 appendMajorVector(vec);
00174 else
00175 appendMinorVector(vec);
00176 }
00177
00178
00179
00180 void
00181 CoinPackedMatrix::appendCol(const int vecsize,
00182 const int *vecind,
00183 const double *vecelem) throw(CoinError)
00184 {
00185 if (colOrdered_)
00186 appendMajorVector(vecsize, vecind, vecelem);
00187 else
00188 appendMinorVector(vecsize, vecind, vecelem);
00189 }
00190
00191
00192
00193 void
00194 CoinPackedMatrix::appendCols(const int numcols,
00195 const CoinPackedVectorBase * const * cols)
00196 throw(CoinError)
00197 {
00198 if (colOrdered_)
00199 appendMajorVectors(numcols, cols);
00200 else
00201 appendMinorVectors(numcols, cols);
00202 }
00203
00204
00205
00206 void
00207 CoinPackedMatrix::appendRow(const CoinPackedVectorBase& vec) throw(CoinError)
00208 {
00209 if (colOrdered_)
00210 appendMinorVector(vec);
00211 else
00212 appendMajorVector(vec);
00213 }
00214
00215
00216
00217 void
00218 CoinPackedMatrix::appendRow(const int vecsize,
00219 const int *vecind,
00220 const double *vecelem) throw(CoinError)
00221 {
00222 if (colOrdered_)
00223 appendMinorVector(vecsize, vecind, vecelem);
00224 else
00225 appendMajorVector(vecsize, vecind, vecelem);
00226 }
00227
00228
00229
00230 void
00231 CoinPackedMatrix::appendRows(const int numrows,
00232 const CoinPackedVectorBase * const * rows)
00233 throw(CoinError)
00234 {
00235 if (colOrdered_)
00236 appendMinorVectors(numrows, rows);
00237 else
00238 appendMajorVectors(numrows, rows);
00239 }
00240
00241
00242
00243 void
00244 CoinPackedMatrix::rightAppendPackedMatrix(const CoinPackedMatrix& matrix)
00245 throw(CoinError)
00246 {
00247 if (colOrdered_) {
00248 if (matrix.colOrdered_) {
00249 majorAppendSameOrdered(matrix);
00250 } else {
00251 majorAppendOrthoOrdered(matrix);
00252 }
00253 } else {
00254 if (matrix.colOrdered_) {
00255 minorAppendOrthoOrdered(matrix);
00256 } else {
00257 minorAppendSameOrdered(matrix);
00258 }
00259 }
00260 }
00261
00262
00263
00264 void
00265 CoinPackedMatrix::bottomAppendPackedMatrix(const CoinPackedMatrix& matrix)
00266 throw(CoinError)
00267 {
00268 if (colOrdered_) {
00269 if (matrix.colOrdered_) {
00270 minorAppendSameOrdered(matrix);
00271 } else {
00272 minorAppendOrthoOrdered(matrix);
00273 }
00274 } else {
00275 if (matrix.colOrdered_) {
00276 majorAppendOrthoOrdered(matrix);
00277 } else {
00278 majorAppendSameOrdered(matrix);
00279 }
00280 }
00281 }
00282
00283
00284
00285 void
00286 CoinPackedMatrix::deleteCols(const int numDel, const int * indDel)
00287 {
00288 if (numDel) {
00289 if (colOrdered_)
00290 deleteMajorVectors(numDel, indDel);
00291 else
00292 deleteMinorVectors(numDel, indDel);
00293 }
00294 }
00295
00296
00297
00298 void
00299 CoinPackedMatrix::deleteRows(const int numDel, const int * indDel)
00300 {
00301 if (numDel) {
00302 if (colOrdered_)
00303 deleteMinorVectors(numDel, indDel);
00304 else
00305 deleteMajorVectors(numDel, indDel);
00306 }
00307 }
00308
00309
00310
00311
00312
00313 void
00314 CoinPackedMatrix::replaceVector(const int index,
00315 const int numReplace,
00316 const double * newElements)
00317 {
00318 if (index >= 0 && index < majorDim_) {
00319 int length = (length_[index] < numReplace) ? length_[index] : numReplace;
00320 CoinMemcpyN(newElements, length, element_ + start_[index]);
00321 } else {
00322 #ifdef COIN_DEBUG
00323 throw CoinError("bad index", "replaceVector", "CoinPackedMatrix");
00324 #endif
00325 }
00326 }
00327
00328
00329
00330 void
00331 CoinPackedMatrix::modifyCoefficient(int row, int column, double newElement,
00332 bool keepZero)
00333 {
00334 int minorIndex,majorIndex;
00335 if (colOrdered_) {
00336 majorIndex=column;
00337 minorIndex=row;
00338 } else {
00339 minorIndex=column;
00340 majorIndex=row;
00341 }
00342 if (majorIndex >= 0 && majorIndex < majorDim_) {
00343 if (minorIndex >= 0 && minorIndex < minorDim_) {
00344 int j;
00345 int end=start_[majorIndex]+length_[majorIndex];;
00346 for (j=start_[majorIndex];j<end;j++) {
00347 if (minorIndex==index_[j]) {
00348
00349 if (newElement||keepZero) {
00350 element_[j]=newElement;
00351 } else {
00352
00353 length_[majorIndex]--;
00354 end--;
00355 size_--;
00356 for (;j<end;j++) {
00357 element_[j]=element_[j+1];
00358 index_[j]=index_[j+1];
00359 }
00360 }
00361 return;
00362 }
00363 }
00364 if (j==end&&(newElement||keepZero)) {
00365
00366 if (end<start_[majorIndex+1]) {
00367 size_++;
00368 length_[majorIndex]++;
00369 element_[end]=newElement;
00370 index_[end]=minorIndex;
00371 } else {
00372 int * addedEntries = new int[majorDim_];
00373 memset(addedEntries, 0, majorDim_ * sizeof(int));
00374 addedEntries[majorIndex] = 1;
00375 resizeForAddingMinorVectors(addedEntries);
00376 delete[] addedEntries;
00377
00378
00379
00380 const int start = start_[majorIndex];
00381 end = start_[majorIndex]+length_[majorIndex];
00382 for (j = end - 1; j >= start; --j) {
00383 if (element_[j] <= newElement)
00384 break;
00385 index_[j+1] = index_[j];
00386 element_[j+1] = element_[j];
00387 }
00388 ++j;
00389 index_[j] = minorIndex;
00390 element_[j] = newElement;
00391 size_++;
00392 length_[majorIndex]++;
00393 }
00394 }
00395 } else {
00396 #ifdef COIN_DEBUG
00397 throw CoinError("bad minor index", "modifyCoefficient",
00398 "CoinPackedMatrix");
00399 #endif
00400 }
00401 } else {
00402 #ifdef COIN_DEBUG
00403 throw CoinError("bad major index", "modifyCoefficient",
00404 "CoinPackedMatrix");
00405 #endif
00406 }
00407 }
00408
00409
00410
00411
00412
00413
00414
00415 int
00416 CoinPackedMatrix::compress(double threshold)
00417 {
00418 CoinBigIndex numberEliminated =0;
00419
00420 int * eliminatedIndex = new int[minorDim_];
00421 double * eliminatedElement = new double[minorDim_];
00422 int i;
00423 for (i=0;i<majorDim_;i++) {
00424 int length = length_[i];
00425 CoinBigIndex k=start_[i];
00426 int kbad=0;
00427 CoinBigIndex j;
00428 for (j=start_[i];j<start_[i]+length;j++) {
00429 if (fabs(element_[j])>=threshold) {
00430 element_[k]=element_[j];
00431 index_[k++]=index_[j];
00432 } else {
00433 eliminatedElement[kbad]=element_[j];
00434 eliminatedIndex[kbad++]=index_[j];
00435 }
00436 }
00437 if (kbad) {
00438 numberEliminated += kbad;
00439 length_[i] = k-start_[i];
00440 memcpy(index_+k,eliminatedIndex,kbad*sizeof(int));
00441 memcpy(element_+k,eliminatedElement,kbad*sizeof(double));
00442 }
00443 }
00444 size_ -= numberEliminated;
00445 delete [] eliminatedIndex;
00446 delete [] eliminatedElement;
00447 return numberEliminated;
00448 }
00449
00450
00451
00452
00453
00454
00455 int
00456 CoinPackedMatrix::eliminateDuplicates(double threshold)
00457 {
00458 CoinBigIndex numberEliminated =0;
00459
00460 int * mark = new int [minorDim_];
00461 int i;
00462 for (i=0;i<minorDim_;i++)
00463 mark[i]=-1;
00464 for (i=0;i<majorDim_;i++) {
00465 CoinBigIndex k=start_[i];
00466 CoinBigIndex end = k+length_[i];
00467 CoinBigIndex j;
00468 for (j=k;j<end;j++) {
00469 int index = index_[j];
00470 if (mark[index]==-1) {
00471 mark[index]=j;
00472 } else {
00473
00474 int jj = mark[index];
00475 element_[jj] += element_[j];
00476 element_[j]=0.0;
00477 }
00478 }
00479 for (j=k;j<end;j++) {
00480 int index = index_[j];
00481 mark[index]=-1;
00482 if (fabs(element_[j])>=threshold) {
00483 element_[k]=element_[j];
00484 index_[k++]=index_[j];
00485 }
00486 }
00487 numberEliminated += end-k;
00488 length_[i] = k-start_[i];
00489 }
00490 size_ -= numberEliminated;
00491 delete [] mark;
00492 return numberEliminated;
00493 }
00494
00495
00496 void
00497 CoinPackedMatrix::removeGaps()
00498 {
00499 for (int i = 1; i < majorDim_; ++i) {
00500 const CoinBigIndex si = start_[i];
00501 const int li = length_[i];
00502 start_[i] = start_[i-1] + length_[i-1];
00503 CoinCopy(index_ + si, index_ + (si + li), index_ + start_[i]);
00504 CoinCopy(element_ + si, element_ + (si + li), element_ + start_[i]);
00505 }
00506 start_[majorDim_] = size_;
00507 }
00508
00509
00510
00511 void
00512 CoinPackedMatrix::submatrixOf(const CoinPackedMatrix& matrix,
00513 const int numMajor, const int * indMajor)
00514 throw(CoinError)
00515 {
00516 int i;
00517
00518 int* sortedIndPtr = CoinTestIndexSet(numMajor, indMajor, matrix.majorDim_,
00519 "submatrixOf");
00520 const int * sortedInd = sortedIndPtr == 0 ? indMajor : sortedIndPtr;
00521
00522 gutsOfDestructor();
00523
00524
00525 CoinBigIndex nzcnt = 0;
00526 const int* length = matrix.getVectorLengths();
00527 for (i = 0; i < numMajor; ++i) {
00528 nzcnt += length[sortedInd[i]];
00529 }
00530
00531 colOrdered_ = matrix.colOrdered_;
00532 maxMajorDim_ = int(numMajor * (1+extraMajor_) + 1);
00533 maxSize_ = int(nzcnt * (1+extraMajor_) * (1+extraGap_) + 100);
00534 length_ = new int[maxMajorDim_];
00535 start_ = new CoinBigIndex[maxMajorDim_+1];
00536 index_ = new int[maxSize_];
00537 element_ = new double[maxSize_];
00538 majorDim_ = 0;
00539 minorDim_ = matrix.minorDim_;
00540 size_ = 0;
00541
00542 for (i = 0; i < numMajor; ++i) {
00543 appendMajorVector(matrix.getVector(sortedInd[i]));
00544 }
00545
00546 delete[] sortedIndPtr;
00547 }
00548
00549
00550
00551 void
00552 CoinPackedMatrix::copyOf(const CoinPackedMatrix& rhs)
00553 {
00554 if (this != &rhs) {
00555 gutsOfDestructor();
00556 gutsOfCopyOf(rhs.colOrdered_,
00557 rhs.minorDim_, rhs.majorDim_, rhs.size_,
00558 rhs.element_, rhs.index_, rhs.start_, rhs.length_,
00559 rhs.extraMajor_, rhs.extraGap_);
00560 }
00561 }
00562
00563
00564
00565 void
00566 CoinPackedMatrix::copyOf(const bool colordered,
00567 const int minor, const int major,
00568 const CoinBigIndex numels,
00569 const double * elem, const int * ind,
00570 const CoinBigIndex * start, const int * len,
00571 const double extraMajor, const double extraGap)
00572 {
00573 gutsOfDestructor();
00574 gutsOfCopyOf(colordered, minor, major, numels, elem, ind, start, len,
00575 extraMajor, extraGap);
00576 }
00577
00578
00579
00580
00581
00582
00583 void
00584 CoinPackedMatrix::reverseOrderedCopyOf(const CoinPackedMatrix& rhs)
00585 {
00586 if (this == &rhs) {
00587 reverseOrdering();
00588 return;
00589 }
00590
00591 int i;
00592 colOrdered_ = !rhs.colOrdered_;
00593 majorDim_ = rhs.minorDim_;
00594 minorDim_ = rhs.majorDim_;
00595 size_ = rhs.size_;
00596
00597 if (size_ == 0) {
00598
00599 maxMajorDim_=majorDim_;
00600 delete[] start_;
00601 delete[] length_;
00602 delete[] index_;
00603 delete[] element_;
00604 start_ = new CoinBigIndex[maxMajorDim_ + 1];
00605 length_ = new int[maxMajorDim_];
00606 for (i = 0; i < majorDim_; ++i) {
00607 start_[i] = 0;
00608 length_[i]=0;
00609 }
00610 start_[majorDim_]=0;
00611 index_ = new int[maxSize_];
00612 element_ = new double[maxSize_];
00613 return;
00614 }
00615
00616
00617
00618
00619 int * orthoLengthPtr = rhs.countOrthoLength();
00620 const int * orthoLength = orthoLengthPtr;
00621
00622
00623
00624 const int newMaxMajorDim_ =
00625 CoinMax(maxMajorDim_, CoinLengthWithExtra(majorDim_, extraMajor_));
00626
00627 if (newMaxMajorDim_ > maxMajorDim_) {
00628 maxMajorDim_ = newMaxMajorDim_;
00629 delete[] start_;
00630 delete[] length_;
00631 start_ = new CoinBigIndex[maxMajorDim_ + 1];
00632 length_ = new int[maxMajorDim_];
00633 }
00634
00635 start_[0] = 0;
00636 if (extraGap_ == 0) {
00637 for (i = 0; i < majorDim_; ++i)
00638 start_[i+1] = start_[i] + orthoLength[i];
00639 } else {
00640 const double eg = extraGap_;
00641 for (i = 0; i < majorDim_; ++i)
00642 start_[i+1] = start_[i] + CoinLengthWithExtra(orthoLength[i], eg);
00643 }
00644
00645 const CoinBigIndex newMaxSize =
00646 CoinMax(maxSize_, getLastStart()+ (CoinBigIndex) extraMajor_);
00647
00648 if (newMaxSize > maxSize_) {
00649 maxSize_ = newMaxSize;
00650 delete[] index_;
00651 delete[] element_;
00652 index_ = new int[maxSize_];
00653 element_ = new double[maxSize_];
00654 }
00655
00656 delete[] orthoLengthPtr;
00657
00658
00659 minorDim_ = 0;
00660 CoinZeroN(length_, majorDim_);
00661 for (i = 0; i < rhs.majorDim_; ++i) {
00662 const CoinBigIndex last = rhs.getVectorLast(i);
00663 for (CoinBigIndex j = rhs.getVectorFirst(i); j != last; ++j) {
00664 const int ind = rhs.index_[j];
00665 element_[start_[ind] + length_[ind]] = rhs.element_[j];
00666 index_[start_[ind] + (length_[ind]++)] = minorDim_;
00667 }
00668 ++minorDim_;
00669 }
00670 }
00671
00672
00673
00674 void
00675 CoinPackedMatrix::assignMatrix(const bool colordered,
00676 const int minor, const int major,
00677 const CoinBigIndex numels,
00678 double *& elem, int *& ind,
00679 CoinBigIndex *& start, int *& len,
00680 const int maxmajor, const CoinBigIndex maxsize)
00681 {
00682 gutsOfDestructor();
00683 colOrdered_ = colordered;
00684 element_ = elem;
00685 index_ = ind;
00686 start_ = start;
00687 majorDim_ = major;
00688 minorDim_ = minor;
00689 size_ = numels;
00690 maxMajorDim_ = maxmajor != -1 ? maxmajor : major;
00691 maxSize_ = maxsize != -1 ? maxsize : numels;
00692 if (len == NULL) {
00693 length_ = new int[maxMajorDim_];
00694 std::adjacent_difference(start + 1, start + (major + 1), length_);
00695 } else {
00696 length_ = len;
00697 }
00698 elem = NULL;
00699 ind = NULL;
00700 start = NULL;
00701 len = NULL;
00702 }
00703
00704
00705
00706 CoinPackedMatrix &
00707 CoinPackedMatrix::operator=(const CoinPackedMatrix& rhs)
00708 {
00709 if (this != &rhs) {
00710 gutsOfDestructor();
00711 gutsOfOpEqual(rhs.colOrdered_,
00712 rhs.minorDim_, rhs.majorDim_, rhs.size_,
00713 rhs.element_, rhs.index_, rhs.start_, rhs.length_);
00714 }
00715 return *this;
00716 }
00717
00718
00719
00720 void
00721 CoinPackedMatrix::reverseOrdering()
00722 {
00723 CoinPackedMatrix m;
00724 m.extraGap_ = extraMajor_;
00725 m.extraMajor_ = extraGap_;
00726 m.reverseOrderedCopyOf(*this);
00727 swap(m);
00728 }
00729
00730
00731
00732 void
00733 CoinPackedMatrix::transpose()
00734 {
00735 colOrdered_ = ! colOrdered_;
00736 }
00737
00738
00739
00740 void
00741 CoinPackedMatrix::swap(CoinPackedMatrix& m)
00742 {
00743 std::swap(colOrdered_, m.colOrdered_);
00744 std::swap(extraGap_, m.extraGap_);
00745 std::swap(extraMajor_, m.extraMajor_);
00746 std::swap(element_, m.element_);
00747 std::swap(index_, m.index_);
00748 std::swap(start_, m.start_);
00749 std::swap(length_, m.length_);
00750 std::swap(majorDim_, m.majorDim_);
00751 std::swap(minorDim_, m.minorDim_);
00752 std::swap(size_, m.size_);
00753 std::swap(maxMajorDim_, m.maxMajorDim_);
00754 std::swap(maxSize_, m.maxSize_);
00755 }
00756
00757
00758
00759
00760 void
00761 CoinPackedMatrix::times(const double * x, double * y) const
00762 {
00763 if (colOrdered_)
00764 timesMajor(x, y);
00765 else
00766 timesMinor(x, y);
00767 }
00768
00769
00770
00771 void
00772 CoinPackedMatrix::times(const CoinPackedVectorBase& x, double * y) const
00773 {
00774 if (colOrdered_)
00775 timesMajor(x, y);
00776 else
00777 timesMinor(x, y);
00778 }
00779
00780
00781
00782 void
00783 CoinPackedMatrix::transposeTimes(const double * x, double * y) const
00784 {
00785 if (colOrdered_)
00786 timesMinor(x, y);
00787 else
00788 timesMajor(x, y);
00789 }
00790
00791
00792
00793 void
00794 CoinPackedMatrix::transposeTimes(const CoinPackedVectorBase& x, double * y) const
00795 {
00796 if (colOrdered_)
00797 timesMinor(x, y);
00798 else
00799 timesMajor(x, y);
00800 }
00801
00802
00803
00804
00805 int *
00806 CoinPackedMatrix::countOrthoLength() const
00807 {
00808 int * orthoLength = new int[minorDim_];
00809 CoinZeroN(orthoLength, minorDim_);
00810 for (int i = majorDim_ - 1; i >= 0; --i) {
00811 const CoinBigIndex last = getVectorLast(i);
00812 for (CoinBigIndex j = getVectorFirst(i); j != last; ++j) {
00813 assert( index_[j] < minorDim_ && index_[j]>=0);
00814 ++orthoLength[index_[j]];
00815 }
00816 }
00817 return orthoLength;
00818 }
00819
00820
00821
00822
00823 void
00824 CoinPackedMatrix::appendMajorVector(const int vecsize,
00825 const int *vecind,
00826 const double *vecelem)
00827 throw(CoinError)
00828 {
00829 #ifdef COIN_DEBUG
00830 for (int i = 0; i < vecsize; ++i) {
00831 if (vecind[i] < 0 || vecind[i] >= minorDim_)
00832 throw CoinError("out of range index",
00833 "appendMajorVector", "CoinPackedMatrix");
00834 }
00835 #if 0
00836 if (std::find_if(vecind, vecind + vecsize,
00837 compose2(logical_or<bool>(),
00838 bind2nd(less<int>(), 0),
00839 bind2nd(greater_equal<int>(), minorDim_))) !=
00840 vecind + vecsize)
00841 throw CoinError("out of range index",
00842 "appendMajorVector", "CoinPackedMatrix");
00843 #endif
00844 #endif
00845
00846 if (majorDim_ == maxMajorDim_ || vecsize > maxSize_ - getLastStart()) {
00847 resizeForAddingMajorVectors(1, &vecsize);
00848 }
00849
00850
00851 const CoinBigIndex last = getLastStart();
00852
00853
00854
00855 length_[majorDim_] = vecsize;
00856 CoinMemcpyN(vecind, vecsize, index_ + last);
00857 CoinMemcpyN(vecelem, vecsize, element_ + last);
00858 if (majorDim_ == 0)
00859 start_[0] = 0;
00860 start_[majorDim_ + 1] =
00861 CoinMin(last + CoinLengthWithExtra(vecsize, extraGap_), maxSize_ );
00862
00863
00864
00865 if (vecsize > 0) {
00866 minorDim_ = CoinMax(minorDim_,
00867 (*std::max_element(vecind, vecind+vecsize)) + 1);
00868 }
00869
00870 ++majorDim_;
00871 size_ += vecsize;
00872 }
00873
00874
00875
00876 void
00877 CoinPackedMatrix::appendMajorVector(const CoinPackedVectorBase& vec)
00878 throw(CoinError)
00879 {
00880 appendMajorVector(vec.getNumElements(),
00881 vec.getIndices(), vec.getElements());
00882 }
00883
00884
00885
00886 void
00887 CoinPackedMatrix::appendMajorVectors(const int numvecs,
00888 const CoinPackedVectorBase * const * vecs)
00889 throw(CoinError)
00890 {
00891 int i;
00892 CoinBigIndex nz = 0;
00893 for (i = 0; i < numvecs; ++i)
00894 nz += CoinLengthWithExtra(vecs[i]->getNumElements(), extraGap_);
00895 reserve(majorDim_ + numvecs, getLastStart() + nz);
00896 for (i = 0; i < numvecs; ++i)
00897 appendMajorVector(*vecs[i]);
00898 }
00899
00900
00901
00902 void
00903 CoinPackedMatrix::appendMinorVector(const int vecsize,
00904 const int *vecind,
00905 const double *vecelem)
00906 throw(CoinError)
00907 {
00908 int i;
00909 #ifdef COIN_DEBUG
00910 for (i = 0; i < vecsize; ++i) {
00911 if (vecind[i] < 0 || vecind[i] >= majorDim_)
00912 throw CoinError("out of range index",
00913 "appendMinorVector", "CoinPackedMatrix");
00914 }
00915 #if 0
00916 if (std::find_if(vecind, vecind + vecsize,
00917 compose2(logical_or<bool>(),
00918 bind2nd(less<int>(), 0),
00919 bind2nd(greater_equal<int>(), majorDim_))) !=
00920 vecind + vecsize)
00921 throw CoinError("out of range index",
00922 "appendMinorVector", "CoinPackedMatrix");
00923 #endif
00924 #endif
00925
00926
00927
00928
00929 for (i = vecsize - 1; i >= 0; --i) {
00930 const int j = vecind[i];
00931 if (start_[j] + length_[j] == start_[j+1])
00932 break;
00933 }
00934
00935 if (i >= 0) {
00936 int * addedEntries = new int[majorDim_];
00937 memset(addedEntries, 0, majorDim_ * sizeof(int));
00938 for (i = vecsize - 1; i >= 0; --i)
00939 addedEntries[vecind[i]] = 1;
00940 resizeForAddingMinorVectors(addedEntries);
00941 delete[] addedEntries;
00942 }
00943
00944
00945 for (i = vecsize - 1; i >= 0; --i) {
00946 const int j = vecind[i];
00947 const CoinBigIndex posj = start_[j] + (length_[j]++);
00948 index_[posj] = minorDim_;
00949 element_[posj] = vecelem[i];
00950 }
00951
00952 ++minorDim_;
00953 size_ += vecsize;
00954 }
00955
00956
00957
00958 void
00959 CoinPackedMatrix::appendMinorVector(const CoinPackedVectorBase& vec)
00960 throw(CoinError)
00961 {
00962 appendMinorVector(vec.getNumElements(),
00963 vec.getIndices(), vec.getElements());
00964 }
00965
00966
00967
00968 void
00969 CoinPackedMatrix::appendMinorVectors(const int numvecs,
00970 const CoinPackedVectorBase * const * vecs)
00971 throw(CoinError)
00972 {
00973 if (numvecs == 0)
00974 return;
00975
00976 int i;
00977
00978 int * addedEntries = new int[majorDim_];
00979 CoinZeroN(addedEntries, majorDim_);
00980 for (i = numvecs - 1; i >= 0; --i) {
00981 const int vecsize = vecs[i]->getNumElements();
00982 const int* vecind = vecs[i]->getIndices();
00983 for (int j = vecsize - 1; j >= 0; --j) {
00984 #ifdef COIN_DEBUG
00985 if (vecind[j] < 0 || vecind[j] >= majorDim_)
00986 throw CoinError("out of range index", "appendMinorVectors",
00987 "CoinPackedMatrix");
00988 #endif
00989 ++addedEntries[vecind[j]];
00990 }
00991 }
00992
00993 for (i = majorDim_ - 1; i >= 0; --i) {
00994 if (start_[i] + length_[i] + addedEntries[i] > start_[i+1])
00995 break;
00996 }
00997 if (i >= 0)
00998 resizeForAddingMinorVectors(addedEntries);
00999 delete[] addedEntries;
01000
01001
01002 for (i = 0; i < numvecs; ++i) {
01003 const int vecsize = vecs[i]->getNumElements();
01004 const int* vecind = vecs[i]->getIndices();
01005 const double* vecelem = vecs[i]->getElements();
01006 for (int j = vecsize - 1; j >= 0; --j) {
01007 const int ind = vecind[j];
01008 element_[start_[ind] + length_[ind]] = vecelem[j];
01009 index_[start_[ind] + (length_[ind]++)] = minorDim_;
01010 }
01011 ++minorDim_;
01012 size_ += vecsize;
01013 }
01014 }
01015
01016
01017
01018
01019 void
01020 CoinPackedMatrix::majorAppendSameOrdered(const CoinPackedMatrix& matrix)
01021 throw(CoinError)
01022 {
01023 if (minorDim_ != matrix.minorDim_) {
01024 throw CoinError("dimension mismatch", "rightAppendSameOrdered",
01025 "CoinPackedMatrix");
01026 }
01027 if (matrix.majorDim_ == 0)
01028 return;
01029
01030 int i;
01031 if (majorDim_ + matrix.majorDim_ > maxMajorDim_ ||
01032 getLastStart() + matrix.getLastStart() > maxSize_) {
01033
01034
01035
01036 resizeForAddingMajorVectors(matrix.majorDim_, matrix.length_);
01037 start_ += majorDim_;
01038 for (i = 0; i < matrix.majorDim_; ++i) {
01039 const int l = matrix.length_[i];
01040 CoinMemcpyN(matrix.index_ + matrix.start_[i], l,
01041 index_ + start_[i]);
01042 CoinMemcpyN(matrix.element_ + matrix.start_[i], l,
01043 element_ + start_[i]);
01044 }
01045 start_ -= majorDim_;
01046 } else {
01047 start_ += majorDim_;
01048 length_ += majorDim_;
01049 for (i = 0; i < matrix.majorDim_; ++i) {
01050 const int l = matrix.length_[i];
01051 CoinMemcpyN(matrix.index_ + matrix.start_[i], l,
01052 index_ + start_[i]);
01053 CoinMemcpyN(matrix.element_ + matrix.start_[i], l,
01054 element_ + start_[i]);
01055 start_[i+1] = start_[i] + matrix.start_[i+1] - matrix.start_[i];
01056 length_[i] = l;
01057 }
01058 start_ -= majorDim_;
01059 length_ -= majorDim_;
01060 }
01061 majorDim_ += matrix.majorDim_;
01062 size_ += matrix.size_;
01063 }
01064
01065
01066
01067 void
01068 CoinPackedMatrix::minorAppendSameOrdered(const CoinPackedMatrix& matrix)
01069 throw(CoinError)
01070 {
01071 if (majorDim_ != matrix.majorDim_) {
01072 throw CoinError("dimension mismatch", "bottomAppendSameOrdered",
01073 "CoinPackedMatrix");
01074 }
01075 if (matrix.minorDim_ == 0)
01076 return;
01077
01078 int i;
01079 for (i = majorDim_ - 1; i >= 0; --i) {
01080 if (start_[i] + length_[i] + matrix.length_[i] > start_[i+1])
01081 break;
01082 }
01083 if (i >= 0)
01084 resizeForAddingMinorVectors(matrix.length_);
01085
01086
01087 for (i = majorDim_ - 1; i >= 0; --i) {
01088 const int l = matrix.length_[i];
01089 std::transform(matrix.index_ + matrix.start_[i],
01090 matrix.index_ + (matrix.start_[i] + l),
01091 index_ + (start_[i] + length_[i]),
01092 std::bind2nd(std::plus<int>(), minorDim_));
01093 CoinMemcpyN(matrix.element_ + matrix.start_[i], l,
01094 element_ + (start_[i] + length_[i]));
01095 length_[i] += l;
01096 }
01097 minorDim_ += matrix.minorDim_;
01098 size_ += matrix.size_;
01099 }
01100
01101
01102
01103 void
01104 CoinPackedMatrix::majorAppendOrthoOrdered(const CoinPackedMatrix& matrix)
01105 throw(CoinError)
01106 {
01107 if (minorDim_ != matrix.majorDim_) {
01108 throw CoinError("dimension mismatch", "rightAppendOrthoOrdered",
01109 "CoinPackedMatrix");
01110 }
01111 if (matrix.majorDim_ == 0)
01112 return;
01113
01114 int i;
01115 CoinBigIndex j;
01116
01117
01118 int * orthoLengthPtr = matrix.countOrthoLength();
01119 const int * orthoLength = orthoLengthPtr;
01120
01121 if (majorDim_ + matrix.minorDim_ > maxMajorDim_) {
01122 resizeForAddingMajorVectors(matrix.minorDim_, orthoLength);
01123 } else {
01124 const double extra_gap = extraGap_;
01125 start_ += majorDim_;
01126 for (i = 0; i < matrix.minorDim_ - 1; ++i) {
01127 start_[i+1] = start_[i] + CoinLengthWithExtra(orthoLength[i], extra_gap);
01128 }
01129 start_ -= majorDim_;
01130 if (start_[majorDim_ + matrix.minorDim_] > maxSize_) {
01131 resizeForAddingMajorVectors(matrix.minorDim_, orthoLength);
01132 }
01133 }
01134
01135
01136
01137
01138
01139
01140 start_ += majorDim_;
01141 length_ += majorDim_;
01142
01143 CoinZeroN(length_, matrix.minorDim_);
01144
01145 for (i = 0; i < matrix.majorDim_; ++i) {
01146 const CoinBigIndex last = matrix.getVectorLast(i);
01147 for (j = matrix.getVectorFirst(i); j < last; ++j) {
01148 const int ind = matrix.index_[j];
01149 element_[start_[ind] + length_[ind]] = matrix.element_[j];
01150 index_[start_[ind] + (length_[ind]++)] = i;
01151 }
01152 }
01153
01154 length_ -= majorDim_;
01155 start_ -= majorDim_;
01156
01157 delete[] orthoLengthPtr;
01158 }
01159
01160
01161
01162 void
01163 CoinPackedMatrix::minorAppendOrthoOrdered(const CoinPackedMatrix& matrix)
01164 throw(CoinError)
01165 {
01166 if (majorDim_ != matrix.minorDim_) {
01167 throw CoinError("dimension mismatch", "bottomAppendOrthoOrdered",
01168 "CoinPackedMatrix");
01169 }
01170 if (matrix.majorDim_ == 0)
01171 return;
01172
01173 int i;
01174
01175
01176
01177
01178 int * addedEntriesPtr = matrix.countOrthoLength();
01179 const int * addedEntries = addedEntriesPtr;
01180 for (i = majorDim_ - 1; i >= 0; --i) {
01181 if (start_[i] + length_[i] + addedEntries[i] > start_[i+1])
01182 break;
01183 }
01184 if (i >= 0)
01185 resizeForAddingMinorVectors(addedEntries);
01186 delete[] addedEntriesPtr;
01187
01188
01189 for (i = 0; i < matrix.majorDim_; ++i) {
01190 const CoinBigIndex last = matrix.getVectorLast(i);
01191 for (CoinBigIndex j = matrix.getVectorFirst(i); j != last; ++j) {
01192 const int ind = matrix.index_[j];
01193 element_[start_[ind] + length_[ind]] = matrix.element_[j];
01194 index_[start_[ind] + (length_[ind]++)] = minorDim_;
01195 }
01196 ++minorDim_;
01197 }
01198 size_ += matrix.size_;
01199 }
01200
01201
01202
01203
01204 void
01205 CoinPackedMatrix::deleteMajorVectors(const int numDel,
01206 const int * indDel) throw(CoinError)
01207 {
01208 int *sortedDelPtr = CoinTestIndexSet(numDel, indDel, majorDim_,
01209 "deleteMajorVectors");
01210 const int * sortedDel = sortedDelPtr == 0 ? indDel : sortedDelPtr;
01211
01212 if (numDel == majorDim_) {
01213
01214 majorDim_ = 0;
01215 minorDim_ = 0;
01216 size_ = 0;
01217 if (sortedDelPtr)
01218 delete[] sortedDelPtr;
01219
01220 maxMajorDim_ = 0;
01221 delete [] length_;
01222 length_ = NULL;
01223 delete [] start_;
01224 start_ = new CoinBigIndex[1];
01225 start_[0]=0;;
01226 delete [] element_;
01227 element_=NULL;
01228 delete [] index_;
01229 index_=NULL;
01230 maxSize_ = 0;
01231 return;
01232 }
01233
01234 CoinBigIndex deleted = 0;
01235 const int last = numDel - 1;
01236 for (int i = 0; i < last; ++i) {
01237 const int ind = sortedDel[i];
01238 const int ind1 = sortedDel[i+1];
01239 deleted += length_[ind];
01240 if (ind1 - ind > 1) {
01241 CoinCopy(start_ + (ind + 1), start_ + ind1, start_ + (ind - i));
01242 CoinCopy(length_ + (ind + 1), length_ + ind1, length_ + (ind - i));
01243 }
01244 }
01245
01246
01247 const int ind = sortedDel[last];
01248 deleted += length_[ind];
01249 if (sortedDel[last] != majorDim_ - 1) {
01250 const int ind1 = majorDim_;
01251 CoinCopy(start_ + (ind + 1), start_ + ind1, start_ + (ind - last));
01252 CoinCopy(length_ + (ind + 1), length_ + ind1, length_ + (ind - last));
01253 }
01254 majorDim_ -= numDel;
01255 const int lastlength = CoinLengthWithExtra(length_[majorDim_-1], extraGap_);
01256 start_[majorDim_] = CoinMin(start_[majorDim_-1] + lastlength, maxSize_);
01257 size_ -= deleted;
01258
01259
01260
01261
01262 if (sortedDel[0] == 0) {
01263 CoinCopyN(index_ + start_[0], length_[0], index_);
01264 CoinCopyN(element_ + start_[0], length_[0], element_);
01265 start_[0] = 0;
01266 }
01267
01268 delete[] sortedDelPtr;
01269 }
01270
01271
01272
01273 void
01274 CoinPackedMatrix::deleteMinorVectors(const int numDel,
01275 const int * indDel) throw(CoinError)
01276 {
01277 if (numDel == minorDim_) {
01278
01279 minorDim_ = 0;
01280 size_ = 0;
01281
01282 memset(length_,0,majorDim_*sizeof(int));
01283 memset(start_,0,(majorDim_+1)*sizeof(int));
01284 delete [] element_;
01285 element_=NULL;
01286 delete [] index_;
01287 index_=NULL;
01288 maxSize_ = 0;
01289 return;
01290 }
01291 int i, j, k;
01292
01293
01294 int* newindexPtr = new int[minorDim_];
01295 CoinIotaN(newindexPtr, minorDim_, 0);
01296 for (j = 0; j < numDel; ++j) {
01297 const int ind = indDel[j];
01298 #ifdef COIN_DEBUG
01299 if (ind < 0 || ind >= minorDim_)
01300 throw CoinError("out of range index",
01301 "deleteMinorVectors", "CoinPackedMatrix");
01302 if (newindexPtr[ind] == -1)
01303 throw CoinError("duplicate index",
01304 "deleteMinorVectors", "CoinPackedMatrix");
01305 #endif
01306 newindexPtr[ind] = -1;
01307 }
01308 for (i = 0, k = 0; i < minorDim_; ++i) {
01309 if (newindexPtr[i] != -1) {
01310 newindexPtr[i] = k++;
01311 }
01312 }
01313
01314 const int * newindex = newindexPtr;
01315 #ifdef TAKEOUT
01316 int mcount[400];
01317 memset(mcount,0,400*sizeof(int));
01318 for (i = 0; i < majorDim_; ++i) {
01319 int * index = index_ + start_[i];
01320 double * elem = element_ + start_[i];
01321 const int length_i = length_[i];
01322 for (j = 0, k = 0; j < length_i; ++j) {
01323 mcount[index[j]]++;
01324 }
01325 }
01326 for (i=0;i<minorDim_;i++) {
01327 if (mcount[i]==10||mcount[i]==15) {
01328 if (newindex[i]>=0)
01329 printf("Keeping original row %d (new %d) with count of %d\n",
01330 i,newindex[i],mcount[i]);
01331 else
01332 printf("deleting row %d with count of %d\n",
01333 i,mcount[i]);
01334 }
01335 }
01336 #endif
01337 int deleted = 0;
01338 for (i = 0; i < majorDim_; ++i) {
01339 int * index = index_ + start_[i];
01340 double * elem = element_ + start_[i];
01341 const int length_i = length_[i];
01342 for (j = 0, k = 0; j < length_i; ++j) {
01343 const int ind = newindex[index[j]];
01344 if (ind != -1) {
01345 index[k] = ind;
01346 elem[k++] = elem[j];
01347 }
01348 }
01349 deleted += length_i - k;
01350 length_[i] = k;
01351 }
01352
01353 delete[] newindexPtr;
01354
01355 minorDim_ -= numDel;
01356 size_ -= deleted;
01357 }
01358
01359
01360
01361
01362 void
01363 CoinPackedMatrix::timesMajor(const double * x, double * y) const
01364 {
01365 memset(y, 0, minorDim_ * sizeof(double));
01366 for (int i = majorDim_ - 1; i >= 0; --i) {
01367 const double x_i = x[i];
01368 if (x_i != 0.0) {
01369 const CoinBigIndex last = getVectorLast(i);
01370 for (CoinBigIndex j = getVectorFirst(i); j < last; ++j)
01371 y[index_[j]] += x_i * element_[j];
01372 }
01373 }
01374 }
01375
01376
01377
01378 void
01379 CoinPackedMatrix::timesMajor(const CoinPackedVectorBase& x, double * y) const
01380 {
01381 memset(y, 0, minorDim_ * sizeof(double));
01382 for (CoinBigIndex i = x.getNumElements() - 1; i >= 0; --i) {
01383 const double x_i = x.getElements()[i];
01384 if (x_i != 0.0) {
01385 const int ind = x.getIndices()[i];
01386 const CoinBigIndex last = getVectorLast(ind);
01387 for (CoinBigIndex j = getVectorFirst(ind); j < last; ++j)
01388 y[index_[j]] += x_i * element_[j];
01389 }
01390 }
01391 }
01392
01393
01394
01395 void
01396 CoinPackedMatrix::timesMinor(const double * x, double * y) const
01397 {
01398 memset(y, 0, majorDim_ * sizeof(double));
01399 for (int i = majorDim_ - 1; i >= 0; --i) {
01400 double y_i = 0;
01401 const CoinBigIndex last = getVectorLast(i);
01402 for (CoinBigIndex j = getVectorFirst(i); j < last; ++j)
01403 y_i += x[index_[j]] * element_[j];
01404 y[i] = y_i;
01405 }
01406 }
01407
01408
01409
01410 void
01411 CoinPackedMatrix::timesMinor(const CoinPackedVectorBase& x, double * y) const
01412 {
01413 memset(y, 0, majorDim_ * sizeof(double));
01414 for (int i = majorDim_ - 1; i >= 0; --i) {
01415 double y_i = 0;
01416 const CoinBigIndex last = getVectorLast(i);
01417 for (CoinBigIndex j = getVectorFirst(i); j < last; ++j)
01418 y_i += x[index_[j]] * element_[j];
01419 y[i] = y_i;
01420 }
01421 }
01422
01423
01424
01425
01426 CoinPackedMatrix::CoinPackedMatrix() :
01427 colOrdered_(true),
01428 extraGap_(0.25),
01429 extraMajor_(0.25),
01430 element_(0),
01431 index_(0),
01432 length_(0),
01433 majorDim_(0),
01434 minorDim_(0),
01435 size_(0),
01436 maxMajorDim_(0),
01437 maxSize_(0)
01438 {
01439 start_ = new CoinBigIndex[1];
01440 start_[0] = 0;
01441 }
01442
01443
01444
01445 CoinPackedMatrix::CoinPackedMatrix(const bool colordered,
01446 const double extraMajor,
01447 const double extraGap) :
01448 colOrdered_(colordered),
01449 extraGap_(extraGap),
01450 extraMajor_(extraMajor),
01451 element_(0),
01452 index_(0),
01453 length_(0),
01454 majorDim_(0),
01455 minorDim_(0),
01456 size_(0),
01457 maxMajorDim_(0),
01458 maxSize_(0)
01459 {
01460 start_ = new CoinBigIndex[1];
01461 start_[0] = 0;
01462 }
01463
01464
01465
01466 CoinPackedMatrix::CoinPackedMatrix(const bool colordered,
01467 const int minor, const int major,
01468 const CoinBigIndex numels,
01469 const double * elem, const int * ind,
01470 const CoinBigIndex * start, const int * len,
01471 const double extraMajor,
01472 const double extraGap) :
01473 colOrdered_(colordered),
01474 extraGap_(extraGap),
01475 extraMajor_(extraMajor),
01476 element_(NULL),
01477 index_(NULL),
01478 start_(NULL),
01479 length_(NULL),
01480 majorDim_(0),
01481 minorDim_(0),
01482 size_(0),
01483 maxMajorDim_(0),
01484 maxSize_(0)
01485 {
01486 gutsOfOpEqual(colordered, minor, major, numels, elem, ind, start, len);
01487 }
01488
01489
01490
01491 CoinPackedMatrix::CoinPackedMatrix(const bool colordered,
01492 const int minor, const int major,
01493 const CoinBigIndex numels,
01494 const double * elem, const int * ind,
01495 const CoinBigIndex * start, const int * len) :
01496 colOrdered_(colordered),
01497 extraGap_(0.0),
01498 extraMajor_(0.0),
01499 element_(NULL),
01500 index_(NULL),
01501 start_(NULL),
01502 length_(NULL),
01503 majorDim_(0),
01504 minorDim_(0),
01505 size_(0),
01506 maxMajorDim_(0),
01507 maxSize_(0)
01508 {
01509 gutsOfOpEqual(colordered, minor, major, numels, elem, ind, start, len);
01510 }
01511
01512
01513
01514
01515
01516
01517
01518
01519
01520
01521
01522
01523
01524
01525
01526
01527
01528
01529
01530
01531
01532
01533
01534
01535
01536
01537
01538
01539
01540
01541
01542
01543
01544
01545
01546
01547
01548 CoinPackedMatrix::CoinPackedMatrix(
01549 const bool colordered,
01550 const int * indexRow ,
01551 const int * indexColumn,
01552 const double * element,
01553 CoinBigIndex numberElements )
01554 :
01555 colOrdered_(colordered),
01556 extraGap_(0.0),
01557 extraMajor_(0.0),
01558 element_(NULL),
01559 index_(NULL),
01560 start_(NULL),
01561 length_(NULL),
01562 majorDim_(0),
01563 minorDim_(0),
01564 size_(0),
01565 maxMajorDim_(0),
01566 maxSize_(0)
01567 {
01568 CoinAbsFltEq eq;
01569 int * colIndices = new int[numberElements];
01570 int * rowIndices = new int[numberElements];
01571 double * elements = new double[numberElements];
01572 CoinCopyN(element,numberElements,elements);
01573 if ( colordered ) {
01574 CoinCopyN(indexColumn,numberElements,colIndices);
01575 CoinCopyN(indexRow,numberElements,rowIndices);
01576 }
01577 else {
01578 CoinCopyN(indexColumn,numberElements,rowIndices);
01579 CoinCopyN(indexRow,numberElements,colIndices);
01580 }
01581
01582 int numberRows=*std::max_element(rowIndices,rowIndices+numberElements)+1;
01583 int * rowCount = new int[numberRows];
01584 int numberColumns=*std::max_element(colIndices,colIndices+numberElements)+1;
01585 int * columnCount = new int[numberColumns];
01586 CoinBigIndex * startColumn = new CoinBigIndex[numberColumns+1];
01587 int * lengths = new int[numberColumns+1];
01588
01589 int iColumn,i;
01590 CoinBigIndex k;
01591 for (i=0;i<numberRows;i++) {
01592 rowCount[i]=0;
01593 }
01594 for (i=0;i<numberColumns;i++) {
01595 columnCount[i]=0;
01596 }
01597 for (i=0;i<numberElements;i++) {
01598 int iRow=rowIndices[i];
01599 int iColumn=colIndices[i];
01600 rowCount[iRow]++;
01601 columnCount[iColumn]++;
01602 }
01603 CoinBigIndex iCount=0;
01604 for (iColumn=0;iColumn<numberColumns;iColumn++) {
01605
01606 iCount+=columnCount[iColumn];
01607 startColumn[iColumn]=iCount;
01608 }
01609 startColumn[iColumn]=iCount;
01610 for (k=numberElements-1;k>=0;k--) {
01611 iColumn=colIndices[k];
01612 if (iColumn>=0) {
01613
01614 double value = elements[k];
01615 int iRow=rowIndices[k];
01616 colIndices[k]=-2;
01617
01618 while (1) {
01619
01620 CoinBigIndex iLook=startColumn[iColumn]-1;
01621 double valueSave=elements[iLook];
01622 int iColumnSave=colIndices[iLook];
01623 int iRowSave=rowIndices[iLook];
01624
01625
01626 startColumn[iColumn]=iLook;
01627 elements[iLook]=value;
01628 rowIndices[iLook]=iRow;
01629 colIndices[iLook]=-1;
01630
01631
01632 if (iColumnSave>=0) {
01633 iColumn=iColumnSave;
01634 value=valueSave;
01635 iRow=iRowSave;
01636 } else if (iColumnSave == -2) {
01637 break;
01638 } else {
01639 assert(1==0);
01640 }
01641
01642 }
01643 }
01644 }
01645
01646
01647
01648 numberElements=0;
01649 for (iColumn=0;iColumn<numberColumns;iColumn++) {
01650 CoinBigIndex start=startColumn[iColumn];
01651 CoinBigIndex end =startColumn[iColumn+1];
01652 lengths[iColumn]=0;
01653 startColumn[iColumn]=numberElements;
01654 if (end>start) {
01655 int lastRow;
01656 double lastValue;
01657
01658 CoinSort_2(rowIndices+start,rowIndices+end,elements+start,CoinFirstLess_2<int, double>());
01659 lastRow=rowIndices[start];
01660 lastValue=elements[start];
01661 for (i=start+1;i<end;i++) {
01662 int iRow=rowIndices[i];
01663 double value=elements[i];
01664 if (iRow>lastRow) {
01665
01666 if(!eq(lastValue,0.0)) {
01667 rowIndices[numberElements]=lastRow;
01668 elements[numberElements]=lastValue;
01669 numberElements++;
01670 lengths[iColumn]++;
01671 }
01672 lastRow=iRow;
01673 lastValue=value;
01674 } else {
01675 lastValue+=value;
01676 }
01677 }
01678
01679 if(!eq(lastValue,0.0)) {
01680 rowIndices[numberElements]=lastRow;
01681 elements[numberElements]=lastValue;
01682 numberElements++;
01683 lengths[iColumn]++;
01684 }
01685 }
01686 }
01687 startColumn[numberColumns]=numberElements;
01688 #if 0
01689 gutsOfOpEqual(colordered,numberRows,numberColumns,numberElements,elements,rowIndices,startColumn,lengths);
01690
01691 delete [] rowCount;
01692 delete [] columnCount;
01693 delete [] startColumn;
01694 delete [] lengths;
01695
01696 delete [] colIndices;
01697 delete [] rowIndices;
01698 delete [] elements;
01699 #else
01700 assignMatrix(colordered,numberRows,numberColumns,numberElements,
01701 elements,rowIndices,startColumn,lengths);
01702 delete [] rowCount;
01703 delete [] columnCount;
01704 delete [] lengths;
01705 delete [] colIndices;
01706 #endif
01707
01708 }
01709
01710
01711
01712 CoinPackedMatrix::CoinPackedMatrix (const CoinPackedMatrix & rhs) :
01713 colOrdered_(true),
01714 extraGap_(0.0),
01715 extraMajor_(0.0),
01716 element_(0),
01717 index_(0),
01718 start_(0),
01719 length_(0),
01720 majorDim_(0),
01721 minorDim_(0),
01722 size_(0),
01723 maxMajorDim_(0),
01724 maxSize_(0)
01725 {
01726 gutsOfCopyOf(rhs.colOrdered_,
01727 rhs.minorDim_, rhs.majorDim_, rhs.size_,
01728 rhs.element_, rhs.index_, rhs.start_, rhs.length_,
01729 rhs.extraMajor_, rhs.extraGap_);
01730 }
01731
01732 CoinPackedMatrix::CoinPackedMatrix (const CoinPackedMatrix & rhs,
01733 int numberRows, const int * whichRow,
01734 int numberColumns,
01735 const int * whichColumn) :
01736 colOrdered_(true),
01737 extraGap_(0.0),
01738 extraMajor_(0.0),
01739 element_(NULL),
01740 index_(NULL),
01741 start_(NULL),
01742 length_(NULL),
01743 majorDim_(0),
01744 minorDim_(0),
01745 size_(0),
01746 maxMajorDim_(0),
01747 maxSize_(0)
01748 {
01749 if (numberRows<=0||numberColumns<=0) {
01750 start_ = new CoinBigIndex[1];
01751 start_[0] = 0;
01752 } else {
01753 if (!rhs.colOrdered_) {
01754
01755 colOrdered_=false;
01756 const int * temp = whichRow;
01757 whichRow = whichColumn;
01758 whichColumn = temp;
01759 int n = numberRows;
01760 numberRows = numberColumns;
01761 numberColumns = n;
01762 }
01763 const double * element1 = rhs.element_;
01764 const int * index1 = rhs.index_;
01765 const CoinBigIndex * start1 = rhs.start_;
01766 const int * length1 = rhs.length_;
01767
01768 majorDim_ = numberColumns;
01769 maxMajorDim_ = numberColumns;
01770 minorDim_ = numberRows;
01771
01772 if (rhs.majorDim_ <= 0 || rhs.minorDim_ <= 0)
01773 throw CoinError("empty rhs", "subset constructor", "CoinPackedMatrix");
01774
01775 int * newRow = new int [rhs.minorDim_];
01776 int iRow;
01777 for (iRow=0;iRow<rhs.minorDim_;iRow++)
01778 newRow[iRow] = -1;
01779
01780 int * duplicateRow = new int [numberRows];
01781 int numberBad=0;
01782 for (iRow=0;iRow<numberRows;iRow++) {
01783 duplicateRow[iRow] = -1;
01784 int kRow = whichRow[iRow];
01785 if (kRow>=0 && kRow < rhs.minorDim_) {
01786 if (newRow[kRow]<0) {
01787
01788 newRow[kRow]=iRow;
01789 } else {
01790
01791 int lastRow = newRow[kRow];
01792 newRow[kRow]=iRow;
01793 duplicateRow[iRow] = lastRow;
01794 }
01795 } else {
01796
01797 numberBad++;
01798 }
01799 }
01800
01801 if (numberBad)
01802 throw CoinError("bad minor entries",
01803 "subset constructor", "CoinPackedMatrix");
01804
01805 size_ = 0;
01806 int iColumn;
01807 numberBad=0;
01808 for (iColumn=0;iColumn<numberColumns;iColumn++) {
01809 int kColumn = whichColumn[iColumn];
01810 if (kColumn>=0 && kColumn <rhs.majorDim_) {
01811 CoinBigIndex i;
01812 for (i=start1[kColumn];i<start1[kColumn]+length1[kColumn];i++) {
01813 int kRow = index1[i];
01814 kRow = newRow[kRow];
01815 while (kRow>=0) {
01816 size_++;
01817 kRow = duplicateRow[kRow];
01818 }
01819 }
01820 } else {
01821
01822 numberBad++;
01823 }
01824 }
01825 if (numberBad)
01826 throw CoinError("bad major entries",
01827 "subset constructor", "CoinPackedMatrix");
01828
01829 maxSize_=max(1,size_);
01830 start_ = new CoinBigIndex [numberColumns+1];
01831 length_ = new int [numberColumns];
01832 index_ = new int[maxSize_];
01833 element_ = new double [maxSize_];
01834
01835 size_ = 0;
01836 start_[0]=0;
01837 for (iColumn=0;iColumn<numberColumns;iColumn++) {
01838 int kColumn = whichColumn[iColumn];
01839 CoinBigIndex i;
01840 for (i=start1[kColumn];i<start1[kColumn]+length1[kColumn];i++) {
01841 int kRow = index1[i];
01842 double value = element1[i];
01843 kRow = newRow[kRow];
01844 while (kRow>=0) {
01845 index_[size_] = kRow;
01846 element_[size_++] = value;
01847 kRow = duplicateRow[kRow];
01848 }
01849 }
01850 start_[iColumn+1] = size_;
01851 length_[iColumn] = size_ - start_[iColumn];
01852 }
01853 delete [] newRow;
01854 delete [] duplicateRow;
01855 }
01856 }
01857
01858
01859
01860
01861 CoinPackedMatrix::~CoinPackedMatrix ()
01862 {
01863 gutsOfDestructor();
01864 }
01865
01866
01867
01868
01869
01870 void
01871 CoinPackedMatrix::gutsOfDestructor()
01872 {
01873 delete[] length_;
01874 delete[] start_;
01875 delete[] index_;
01876 delete[] element_;
01877 length_ = 0;
01878 start_ = 0;
01879 index_ = 0;
01880 element_ = 0;
01881 }
01882
01883
01884
01885 void
01886 CoinPackedMatrix::gutsOfCopyOf(const bool colordered,
01887 const int minor, const int major,
01888 const CoinBigIndex numels,
01889 const double * elem, const int * ind,
01890 const CoinBigIndex * start, const int * len,
01891 const double extraMajor, const double extraGap)
01892 {
01893 colOrdered_ = colordered;
01894 majorDim_ = major;
01895 minorDim_ = minor;
01896 size_ = numels;
01897
01898 extraGap_ = extraGap;
01899 extraMajor_ = extraMajor;
01900
01901 maxMajorDim_ = CoinLengthWithExtra(majorDim_, extraMajor_);
01902
01903 if (maxMajorDim_ > 0) {
01904 length_ = new int[maxMajorDim_];
01905 if (len == 0) {
01906 std::adjacent_difference(start + 1, start + (major + 1), length_);
01907 } else {
01908 CoinMemcpyN(len, major, length_);
01909 }
01910 start_ = new CoinBigIndex[maxMajorDim_+1];
01911 CoinMemcpyN(start, major+1, start_);
01912 }
01913
01914 maxSize_ = maxMajorDim_ > 0 ? start_[major] : 0;
01915 maxSize_ = CoinLengthWithExtra(maxSize_, extraMajor_);
01916
01917 if (maxSize_ > 0) {
01918 element_ = new double[maxSize_];
01919 index_ = new int[maxSize_];
01920
01921
01922
01923 for (int i = majorDim_ - 1; i >= 0; --i) {
01924 CoinMemcpyN(ind + start[i], length_[i], index_ + start_[i]);
01925 CoinMemcpyN(elem + start[i], length_[i], element_ + start_[i]);
01926 }
01927 }
01928 }
01929
01930
01931
01932 void
01933 CoinPackedMatrix::gutsOfOpEqual(const bool colordered,
01934 const int minor, const int major,
01935 const CoinBigIndex numels,
01936 const double * elem, const int * ind,
01937 const CoinBigIndex * start, const int * len)
01938 {
01939 colOrdered_ = colordered;
01940 majorDim_ = major;
01941 minorDim_ = minor;
01942 size_ = numels;
01943
01944 maxMajorDim_ = CoinLengthWithExtra(majorDim_, extraMajor_);
01945
01946 int i;
01947 if (maxMajorDim_ > 0) {
01948 length_ = new int[maxMajorDim_];
01949 if (len == 0) {
01950 std::adjacent_difference(start + 1, start + (major + 1), length_);
01951 } else {
01952 CoinMemcpyN(len, major, length_);
01953 }
01954 start_ = new CoinBigIndex[maxMajorDim_+1];
01955 start_[0] = 0;
01956 if (extraGap_ == 0) {
01957 for (i = 0; i < major; ++i)
01958 start_[i+1] = start_[i] + length_[i];
01959 } else {
01960 const double extra_gap = extraGap_;
01961 for (i = 0; i < major; ++i)
01962 start_[i+1] = start_[i] + CoinLengthWithExtra(length_[i], extra_gap);
01963 }
01964 } else {
01965
01966 start_ = new CoinBigIndex[1];
01967 start_[0] = 0;
01968 }
01969
01970 maxSize_ = maxMajorDim_ > 0 ? start_[major] : 0;
01971 maxSize_ = CoinLengthWithExtra(maxSize_, extraMajor_);
01972
01973 if (maxSize_ > 0) {
01974 element_ = new double[maxSize_];
01975 index_ = new int[maxSize_];
01976
01977
01978
01979 for (i = majorDim_ - 1; i >= 0; --i) {
01980 CoinMemcpyN(ind + start[i], length_[i], index_ + start_[i]);
01981 CoinMemcpyN(elem + start[i], length_[i], element_ + start_[i]);
01982 }
01983 }
01984 }
01985
01986
01987
01988
01989 void
01990 CoinPackedMatrix::resizeForAddingMajorVectors(const int numVec,
01991 const int * lengthVec)
01992 {
01993 const double extra_gap = extraGap_;
01994 int i;
01995
01996 maxMajorDim_ =
01997 CoinMax(maxMajorDim_, CoinLengthWithExtra(majorDim_ + numVec, extraMajor_));
01998
01999 CoinBigIndex * newStart = new CoinBigIndex[maxMajorDim_ + 1];
02000 int * newLength = new int[maxMajorDim_];
02001
02002 CoinMemcpyN(length_, majorDim_, newLength);
02003
02004 CoinMemcpyN(lengthVec, numVec, newLength + majorDim_);
02005 majorDim_ += numVec;
02006
02007 newStart[0] = 0;
02008 if (extra_gap == 0) {
02009 for (i = 0; i < majorDim_; ++i)
02010 newStart[i+1] = newStart[i] + newLength[i];
02011 } else {
02012 for (i = 0; i < majorDim_; ++i)
02013 newStart[i+1] = newStart[i] + CoinLengthWithExtra(newLength[i],extra_gap);
02014 }
02015
02016 maxSize_ =
02017 CoinMax(maxSize_,newStart[majorDim_]+ (int) extraMajor_);
02018 majorDim_ -= numVec;
02019
02020 int * newIndex = new int[maxSize_];
02021 double * newElem = new double[maxSize_];
02022 for (i = majorDim_ - 1; i >= 0; --i) {
02023 CoinMemcpyN(index_ + start_[i], length_[i], newIndex + newStart[i]);
02024 CoinMemcpyN(element_ + start_[i], length_[i], newElem + newStart[i]);
02025 }
02026
02027 gutsOfDestructor();
02028 start_ = newStart;
02029 length_ = newLength;
02030 index_ = newIndex;
02031 element_ = newElem;
02032 }
02033
02034
02035
02036 void
02037 CoinPackedMatrix::resizeForAddingMinorVectors(const int * addedEntries)
02038 {
02039 int i;
02040 maxMajorDim_ =
02041 CoinMax(CoinLengthWithExtra(majorDim_, extraMajor_), maxMajorDim_);
02042 CoinBigIndex * newStart = new CoinBigIndex[maxMajorDim_ + 1];
02043 int * newLength = new int[maxMajorDim_];
02044
02045
02046
02047 for (i = majorDim_ - 1; i >= 0; --i)
02048 newLength[i] = length_[i] + addedEntries[i];
02049
02050 newStart[0] = 0;
02051 if (extraGap_ == 0) {
02052 for (i = 0; i < majorDim_; ++i)
02053 newStart[i+1] = newStart[i] + newLength[i];
02054 } else {
02055 const double eg = extraGap_;
02056 for (i = 0; i < majorDim_; ++i)
02057 newStart[i+1] = newStart[i] + CoinLengthWithExtra(newLength[i], eg);
02058 }
02059
02060
02061 for (i = majorDim_ - 1; i >= 0; --i)
02062 newLength[i] -= addedEntries[i];
02063
02064 maxSize_ =
02065 CoinMax(newStart[majorDim_]+ (int) extraMajor_, maxSize_);
02066 int * newIndex = new int[maxSize_];
02067 double * newElem = new double[maxSize_];
02068 for (i = majorDim_ - 1; i >= 0; --i) {
02069 CoinMemcpyN(index_ + start_[i], length_[i],
02070 newIndex + newStart[i]);
02071 CoinMemcpyN(element_ + start_[i], length_[i],
02072 newElem + newStart[i]);
02073 }
02074
02075 gutsOfDestructor();
02076 start_ = newStart;
02077 length_ = newLength;
02078 index_ = newIndex;
02079 element_ = newElem;
02080 }
02081
02082
02083
02084
02085 void
02086 CoinPackedMatrix::dumpMatrix(const char* fname) const
02087 {
02088 if (! fname) {
02089 printf("Dumping matrix...\n\n");
02090 printf("colordered: %i\n", isColOrdered() ? 1 : 0);
02091 const int major = getMajorDim();
02092 const int minor = getMinorDim();
02093 printf("major: %i minor: %i\n", major, minor);
02094 for (int i = 0; i < major; ++i) {
02095 printf("vec %i has length %i with entries:\n", i, length_[i]);
02096 for (CoinBigIndex j = start_[i]; j < start_[i] + length_[i]; ++j) {
02097 printf(" %15i %40.25f\n", index_[j], element_[j]);
02098 }
02099 }
02100 printf("\nFinished dumping matrix\n");
02101 } else {
02102 FILE* out = fopen(fname, "w");
02103 fprintf(out, "Dumping matrix...\n\n");
02104 fprintf(out, "colordered: %i\n", isColOrdered() ? 1 : 0);
02105 const int major = getMajorDim();
02106 const int minor = getMinorDim();
02107 fprintf(out, "major: %i minor: %i\n", major, minor);
02108 for (int i = 0; i < major; ++i) {
02109 fprintf(out, "vec %i has length %i with entries:\n", i, length_[i]);
02110 for (CoinBigIndex j = start_[i]; j < start_[i] + length_[i]; ++j) {
02111 fprintf(out, " %15i %40.25f\n", index_[j], element_[j]);
02112 }
02113 }
02114 fprintf(out, "\nFinished dumping matrix\n");
02115 fclose(out);
02116 }
02117 }
02118 bool
02119 CoinPackedMatrix::isEquivalent2(const CoinPackedMatrix& rhs) const
02120 {
02121 CoinRelFltEq eq;
02122
02123 if (isColOrdered() ^ rhs.isColOrdered()) {
02124 std::cerr<<"Ordering "<<isColOrdered()<<
02125 " rhs - "<<rhs.isColOrdered()<<std::endl;
02126 return false;
02127 }
02128 if (getNumCols() != rhs.getNumCols()) {
02129 std::cerr<<"NumCols "<<getNumCols()<<
02130 " rhs - "<<rhs.getNumCols()<<std::endl;
02131 return false;
02132 }
02133 if (getNumRows() != rhs.getNumRows()) {
02134 std::cerr<<"NumRows "<<getNumRows()<<
02135 " rhs - "<<rhs.getNumRows()<<std::endl;
02136 return false;
02137 }
02138 if (getNumElements() != rhs.getNumElements()) {
02139 std::cerr<<"NumElements "<<getNumElements()<<
02140 " rhs - "<<rhs.getNumElements()<<std::endl;
02141 return false;
02142 }
02143
02144 for (int i=getMajorDim()-1; i >= 0; --i) {
02145 CoinShallowPackedVector pv = getVector(i);
02146 CoinShallowPackedVector rhsPv = rhs.getVector(i);
02147 if ( !pv.isEquivalent(rhsPv,eq) ) {
02148 std::cerr<<"vector # "<<i<<" nel "<<pv.getNumElements()<<
02149 " rhs - "<<rhsPv.getNumElements()<<std::endl;
02150 int j;
02151 const int * inds = pv.getIndices();
02152 const double * elems = pv.getElements();
02153 const int * inds2 = rhsPv.getIndices();
02154 const double * elems2 = rhsPv.getElements();
02155 for ( j = 0 ;j < pv.getNumElements() ; ++j) {
02156 double diff = elems[j]-elems2[j];
02157 if (diff) {
02158 std::cerr<<j<<"( "<<inds[j]<<", "<<elems[j]<<"), rhs ( "<<
02159 inds2[j]<<", "<<elems2[j]<<") diff "<<
02160 diff<<std::endl;
02161 const int * xx = (const int *) (elems+j);
02162 printf("%x %x",xx[0],xx[1]);
02163 xx = (const int *) (elems2+j);
02164 printf(" %x %x\n",xx[0],xx[1]);
02165 }
02166 }
02167
02168 }
02169 }
02170 return true;
02171 }