00001
00002
00003
00004
00005 #include <cstdio>
00006
00007 #include "CoinPragma.hpp"
00008 #include "CoinIndexedVector.hpp"
00009 #include "CoinHelperFunctions.hpp"
00010
00011 #include "ClpSimplex.hpp"
00012 #include "ClpFactorization.hpp"
00013
00014 #include "ClpPlusMinusOneMatrix.hpp"
00015 #include "ClpMessage.hpp"
00016
00017
00018
00019
00020
00021
00022
00023
00024 ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix ()
00025 : ClpMatrixBase()
00026 {
00027 setType(12);
00028 elements_ = NULL;
00029 startPositive_ = NULL;
00030 startNegative_ = NULL;
00031 lengths_=NULL;
00032 indices_=NULL;
00033 numberRows_=0;
00034 numberColumns_=0;
00035 columnOrdered_=true;
00036 }
00037
00038
00039
00040
00041 ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const ClpPlusMinusOneMatrix & rhs)
00042 : ClpMatrixBase(rhs)
00043 {
00044 elements_ = NULL;
00045 startPositive_ = NULL;
00046 startNegative_ = NULL;
00047 lengths_=NULL;
00048 indices_=NULL;
00049 numberRows_=rhs.numberRows_;
00050 numberColumns_=rhs.numberColumns_;
00051 columnOrdered_=rhs.columnOrdered_;
00052 if (numberColumns_) {
00053 int numberElements = rhs.startPositive_[numberColumns_];
00054 indices_ = new int [ numberElements];
00055 memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
00056 startPositive_ = new int [ numberColumns_+1];
00057 memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
00058 startNegative_ = new int [ numberColumns_];
00059 memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(int));
00060 }
00061 }
00062
00063 ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const CoinPackedMatrix & rhs)
00064 : ClpMatrixBase()
00065 {
00066 setType(12);
00067 elements_ = NULL;
00068 startPositive_ = NULL;
00069 startNegative_ = NULL;
00070 lengths_=NULL;
00071 indices_=NULL;
00072 int iColumn;
00073 assert (rhs.isColOrdered());
00074
00075 const int * row = rhs.getIndices();
00076 const CoinBigIndex * columnStart = rhs.getVectorStarts();
00077 const int * columnLength = rhs.getVectorLengths();
00078 const double * elementByColumn = rhs.getElements();
00079 numberColumns_ = rhs.getNumCols();
00080 bool goodPlusMinusOne=true;
00081 numberRows_=-1;
00082 indices_ = new int[rhs.getNumElements()];
00083 startPositive_ = new int [numberColumns_+1];
00084 startNegative_ = new int [numberColumns_];
00085 int * temp = new int [rhs.getNumRows()];
00086 CoinBigIndex j=0;
00087 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
00088 CoinBigIndex k;
00089 int iNeg=0;
00090 startPositive_[iColumn]=j;
00091 for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];
00092 k++) {
00093 int iRow;
00094 if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
00095 iRow = row[k];
00096 numberRows_ = max(numberRows_,iRow);
00097 indices_[j++]=iRow;
00098 } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
00099 iRow = row[k];
00100 numberRows_ = max(numberRows_,iRow);
00101 temp[iNeg++]=iRow;
00102 } else {
00103 goodPlusMinusOne = false;
00104 }
00105 }
00106 if (goodPlusMinusOne) {
00107
00108 startNegative_[iColumn]=j;
00109 for (k=0;k<iNeg;k++) {
00110 indices_[j++] = temp[k];
00111 }
00112 } else {
00113 break;
00114 }
00115 }
00116 startPositive_[numberColumns_]=j;
00117 delete [] temp;
00118 if (!goodPlusMinusOne) {
00119 delete [] indices_;
00120
00121 indices_=NULL;
00122 numberRows_=0;
00123 numberColumns_=0;
00124 delete [] startPositive_;
00125 delete [] startNegative_;
00126 startPositive_ = NULL;
00127 startNegative_ = NULL;
00128 } else {
00129 numberRows_ ++;
00130 columnOrdered_ = true;
00131 }
00132 }
00133
00134
00135
00136
00137 ClpPlusMinusOneMatrix::~ClpPlusMinusOneMatrix ()
00138 {
00139 delete [] elements_;
00140 delete [] startPositive_;
00141 delete [] startNegative_;
00142 delete [] lengths_;
00143 delete [] indices_;
00144 }
00145
00146
00147
00148
00149 ClpPlusMinusOneMatrix &
00150 ClpPlusMinusOneMatrix::operator=(const ClpPlusMinusOneMatrix& rhs)
00151 {
00152 if (this != &rhs) {
00153 ClpMatrixBase::operator=(rhs);
00154 delete [] elements_;
00155 delete [] startPositive_;
00156 delete [] startNegative_;
00157 delete [] lengths_;
00158 delete [] indices_;
00159 elements_ = NULL;
00160 startPositive_ = NULL;
00161 lengths_=NULL;
00162 indices_=NULL;
00163 numberRows_=rhs.numberRows_;
00164 numberColumns_=rhs.numberColumns_;
00165 columnOrdered_=rhs.columnOrdered_;
00166 if (numberColumns_) {
00167 int numberElements = rhs.startPositive_[numberColumns_];
00168 indices_ = new int [ numberElements];
00169 memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
00170 startPositive_ = new int [ numberColumns_+1];
00171 memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
00172 startNegative_ = new int [ numberColumns_];
00173 memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(int));
00174 }
00175 }
00176 return *this;
00177 }
00178
00179
00180
00181 ClpMatrixBase * ClpPlusMinusOneMatrix::clone() const
00182 {
00183 return new ClpPlusMinusOneMatrix(*this);
00184 }
00185
00186
00187 ClpMatrixBase *
00188 ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
00189 int numberColumns,
00190 const int * whichColumns) const
00191 {
00192 return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
00193 numberColumns, whichColumns);
00194 }
00195
00196
00197 ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
00198 const ClpPlusMinusOneMatrix & rhs,
00199 int numberRows, const int * whichRow,
00200 int numberColumns, const int * whichColumn)
00201 : ClpMatrixBase(rhs)
00202 {
00203 elements_ = NULL;
00204 startPositive_ = NULL;
00205 startNegative_ = NULL;
00206 lengths_=NULL;
00207 indices_=NULL;
00208 numberRows_=0;
00209 numberColumns_=0;
00210 columnOrdered_=rhs.columnOrdered_;
00211 if (numberRows<=0||numberColumns<=0) {
00212 startPositive_ = new int[1];
00213 startPositive_[0] = 0;
00214 } else {
00215 numberColumns_ = numberColumns;
00216 numberRows_ = numberRows;
00217 const int * index1 = rhs.indices_;
00218 int * startPositive1 = rhs.startPositive_;
00219
00220 int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
00221 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
00222 int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
00223 int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
00224
00225 if (!columnOrdered_) {
00226 int temp1 = numberRows;
00227 numberRows = numberColumns;
00228 numberColumns = temp1;
00229 const int * temp2;
00230 temp2 = whichRow;
00231 whichRow = whichColumn;
00232 whichColumn = temp2;
00233 }
00234
00235 if (numberMajor1 <= 0 || numberMinor1 <= 0)
00236 throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
00237
00238 int * newRow = new int [numberMinor1];
00239 int iRow;
00240 for (iRow=0;iRow<numberMinor1;iRow++)
00241 newRow[iRow] = -1;
00242
00243 int * duplicateRow = new int [numberMinor];
00244 int numberBad=0;
00245 for (iRow=0;iRow<numberMinor;iRow++) {
00246 duplicateRow[iRow] = -1;
00247 int kRow = whichRow[iRow];
00248 if (kRow>=0 && kRow < numberMinor1) {
00249 if (newRow[kRow]<0) {
00250
00251 newRow[kRow]=iRow;
00252 } else {
00253
00254 int lastRow = newRow[kRow];
00255 newRow[kRow]=iRow;
00256 duplicateRow[iRow] = lastRow;
00257 }
00258 } else {
00259
00260 numberBad++;
00261 }
00262 }
00263
00264 if (numberBad)
00265 throw CoinError("bad minor entries",
00266 "subset constructor", "ClpPlusMinusOneMatrix");
00267
00268 int size = 0;
00269 int iColumn;
00270 numberBad=0;
00271 for (iColumn=0;iColumn<numberMajor;iColumn++) {
00272 int kColumn = whichColumn[iColumn];
00273 if (kColumn>=0 && kColumn <numberMajor1) {
00274 int i;
00275 for (i=startPositive1[kColumn];i<startPositive1[kColumn+1];i++) {
00276 int kRow = index1[i];
00277 kRow = newRow[kRow];
00278 while (kRow>=0) {
00279 size++;
00280 kRow = duplicateRow[kRow];
00281 }
00282 }
00283 } else {
00284
00285 numberBad++;
00286 printf("%d %d %d %d\n",iColumn,numberMajor,numberMajor1,kColumn);
00287 }
00288 }
00289 if (numberBad)
00290 throw CoinError("bad major entries",
00291 "subset constructor", "ClpPlusMinusOneMatrix");
00292
00293 startPositive_ = new int [numberMajor+1];
00294 startNegative_ = new int [numberMajor];
00295 indices_ = new int[size];
00296
00297 size = 0;
00298 startPositive_[0]=0;
00299 int * startNegative1 = rhs.startNegative_;
00300 for (iColumn=0;iColumn<numberMajor;iColumn++) {
00301 int kColumn = whichColumn[iColumn];
00302 int i;
00303 for (i=startPositive1[kColumn];i<startNegative1[kColumn];i++) {
00304 int kRow = index1[i];
00305 kRow = newRow[kRow];
00306 while (kRow>=0) {
00307 indices_[size++] = kRow;
00308 kRow = duplicateRow[kRow];
00309 }
00310 }
00311 startNegative_[iColumn] = size;
00312 for (;i<startPositive1[kColumn+1];i++) {
00313 int kRow = index1[i];
00314 kRow = newRow[kRow];
00315 while (kRow>=0) {
00316 indices_[size++] = kRow;
00317 kRow = duplicateRow[kRow];
00318 }
00319 }
00320 startPositive_[iColumn+1] = size;
00321 }
00322 delete [] newRow;
00323 delete [] duplicateRow;
00324 }
00325 }
00326
00327
00328
00329 ClpMatrixBase *
00330 ClpPlusMinusOneMatrix::reverseOrderedCopy() const
00331 {
00332 int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
00333 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
00334
00335 int * tempP = new int [numberMinor];
00336 int * tempN = new int [numberMinor];
00337 memset(tempP,0,numberMinor*sizeof(int));
00338 memset(tempN,0,numberMinor*sizeof(int));
00339 CoinBigIndex j=0;
00340 int i;
00341 for (i=0;i<numberMajor;i++) {
00342 for (;j<startNegative_[i];j++) {
00343 int iRow = indices_[j];
00344 tempP[iRow]++;
00345 }
00346 for (;j<startPositive_[i+1];j++) {
00347 int iRow = indices_[j];
00348 tempN[iRow]++;
00349 }
00350 }
00351 int * newIndices = new int [startPositive_[numberMajor]];
00352 int * newP = new int [numberMinor+1];
00353 int * newN = new int[numberMinor];
00354 int iRow;
00355 j=0;
00356
00357 for (iRow=0;iRow<numberMinor;iRow++) {
00358 newP[iRow]=j;
00359 j += tempP[iRow];
00360 tempP[iRow]=newP[iRow];
00361 newN[iRow] = j;
00362 j += tempN[iRow];
00363 tempN[iRow]=newN[iRow];
00364 }
00365 newP[numberMinor]=j;
00366 j=0;
00367 for (i=0;i<numberMajor;i++) {
00368 for (;j<startNegative_[i];j++) {
00369 int iRow = indices_[j];
00370 int put = tempP[iRow];
00371 newIndices[put++] = i;
00372 tempP[iRow] = put;
00373 }
00374 for (;j<startPositive_[i+1];j++) {
00375 int iRow = indices_[j];
00376 int put = tempN[iRow];
00377 newIndices[put++] = i;
00378 tempN[iRow] = put;
00379 }
00380 }
00381 delete [] tempP;
00382 delete [] tempN;
00383 ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
00384 newCopy->passInCopy(numberMinor, numberMajor,
00385 !columnOrdered_, newIndices, newP, newN);
00386 return newCopy;
00387 }
00388
00389 void
00390 ClpPlusMinusOneMatrix::times(double scalar,
00391 const double * x, double * y) const
00392 {
00393 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
00394 int i;
00395 CoinBigIndex j;
00396 assert (columnOrdered_);
00397 for (i=0;i<numberMajor;i++) {
00398 double value = scalar*x[i];
00399 if (value) {
00400 for (j=startPositive_[i];j<startNegative_[i];j++) {
00401 int iRow = indices_[j];
00402 y[iRow] += value;
00403 }
00404 for (;j<startPositive_[i+1];j++) {
00405 int iRow = indices_[j];
00406 y[iRow] -= value;
00407 }
00408 }
00409 }
00410 }
00411 void
00412 ClpPlusMinusOneMatrix::transposeTimes(double scalar,
00413 const double * x, double * y) const
00414 {
00415 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
00416 int i;
00417 CoinBigIndex j=0;
00418 assert (columnOrdered_);
00419 for (i=0;i<numberMajor;i++) {
00420 double value = 0.0;
00421 for (;j<startNegative_[i];j++) {
00422 int iRow = indices_[j];
00423 value += x[iRow];
00424 }
00425 for (;j<startPositive_[i+1];j++) {
00426 int iRow = indices_[j];
00427 value -= x[iRow];
00428 }
00429 y[i] += scalar*value;
00430 }
00431 }
00432 void
00433 ClpPlusMinusOneMatrix::times(double scalar,
00434 const double * x, double * y,
00435 const double * rowScale,
00436 const double * columnScale) const
00437 {
00438
00439 times(scalar, x, y);
00440 }
00441 void
00442 ClpPlusMinusOneMatrix::transposeTimes( double scalar,
00443 const double * x, double * y,
00444 const double * rowScale,
00445 const double * columnScale) const
00446 {
00447
00448 transposeTimes(scalar, x, y);
00449 }
00450
00451
00452 void
00453 ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex * model, double scalar,
00454 const CoinIndexedVector * rowArray,
00455 CoinIndexedVector * y,
00456 CoinIndexedVector * columnArray) const
00457 {
00458
00459 columnArray->clear();
00460 double * pi = rowArray->denseVector();
00461 int numberNonZero=0;
00462 int * index = columnArray->getIndices();
00463 double * array = columnArray->denseVector();
00464 int numberInRowArray = rowArray->getNumElements();
00465
00466 double zeroTolerance = model->factorization()->zeroTolerance();
00467 int numberRows = model->numberRows();
00468 bool packed = rowArray->packedMode();
00469 ClpPlusMinusOneMatrix* rowCopy =
00470 dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
00471 double factor = 0.3;
00472
00473 int numberColumns = model->numberColumns();
00474
00475
00476 if (numberColumns*sizeof(double)>1000000) {
00477 if (numberRows*10<numberColumns)
00478 factor=0.1;
00479 else if (numberRows*4<numberColumns)
00480 factor=0.15;
00481 else if (numberRows*2<numberColumns)
00482 factor=0.2;
00483 }
00484 if (numberInRowArray>factor*numberRows||!rowCopy) {
00485 assert (!y->getNumElements());
00486
00487
00488 int iColumn;
00489 CoinBigIndex j=0;
00490 assert (columnOrdered_);
00491 if (packed) {
00492
00493 assert(y->capacity()>=numberRows);
00494 double * piOld = pi;
00495 pi = y->denseVector();
00496 const int * whichRow = rowArray->getIndices();
00497 int i;
00498
00499 for (i=0;i<numberInRowArray;i++) {
00500 int iRow = whichRow[i];
00501 pi[iRow]=scalar*piOld[i];
00502 }
00503 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
00504 double value = 0.0;
00505 for (;j<startNegative_[iColumn];j++) {
00506 int iRow = indices_[j];
00507 value += pi[iRow];
00508 }
00509 for (;j<startPositive_[iColumn+1];j++) {
00510 int iRow = indices_[j];
00511 value -= pi[iRow];
00512 }
00513 if (fabs(value)>zeroTolerance) {
00514 array[numberNonZero]=value;
00515 index[numberNonZero++]=iColumn;
00516 }
00517 }
00518 for (i=0;i<numberInRowArray;i++) {
00519 int iRow = whichRow[i];
00520 pi[iRow]=0.0;
00521 }
00522 } else {
00523 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
00524 double value = 0.0;
00525 for (;j<startNegative_[iColumn];j++) {
00526 int iRow = indices_[j];
00527 value += pi[iRow];
00528 }
00529 for (;j<startPositive_[iColumn+1];j++) {
00530 int iRow = indices_[j];
00531 value -= pi[iRow];
00532 }
00533 value *= scalar;
00534 if (fabs(value)>zeroTolerance) {
00535 index[numberNonZero++]=iColumn;
00536 array[iColumn]=value;
00537 }
00538 }
00539 }
00540 columnArray->setNumElements(numberNonZero);
00541 } else {
00542
00543 rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
00544 }
00545 }
00546
00547
00548 void
00549 ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
00550 const CoinIndexedVector * rowArray,
00551 CoinIndexedVector * y,
00552 CoinIndexedVector * columnArray) const
00553 {
00554 columnArray->clear();
00555 double * pi = rowArray->denseVector();
00556 int numberNonZero=0;
00557 int * index = columnArray->getIndices();
00558 double * array = columnArray->denseVector();
00559 int numberInRowArray = rowArray->getNumElements();
00560
00561 double zeroTolerance = model->factorization()->zeroTolerance();
00562 const int * column = indices_;
00563 const CoinBigIndex * startPositive = startPositive_;
00564 const CoinBigIndex * startNegative = startNegative_;
00565 const int * whichRow = rowArray->getIndices();
00566 bool packed = rowArray->packedMode();
00567 if (numberInRowArray>2||y->getNumElements()) {
00568
00569 int iRow;
00570 double * markVector = y->denseVector();
00571 int * mark = y->getIndices();
00572 int numberOriginal=y->getNumElements();
00573 int i;
00574 if (packed) {
00575 assert(!numberOriginal);
00576 numberNonZero=0;
00577
00578 char * marked = (char *) (index+columnArray->capacity());
00579 double * array2 = y->denseVector();
00580 #ifdef CLP_DEBUG
00581 int numberColumns = model->numberColumns();
00582 for (i=0;i<numberColumns;i++) {
00583 assert(!marked[i]);
00584 assert(!array2[i]);
00585 }
00586 #endif
00587 for (i=0;i<numberInRowArray;i++) {
00588 iRow = whichRow[i];
00589 double value = pi[i]*scalar;
00590 CoinBigIndex j;
00591 for (j=startPositive[iRow];j<startNegative[iRow];j++) {
00592 int iColumn = column[j];
00593 if (!marked[iColumn]) {
00594 marked[iColumn]=1;
00595 index[numberNonZero++]=iColumn;
00596 }
00597 array2[iColumn] += value;
00598 }
00599 for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
00600 int iColumn = column[j];
00601 if (!marked[iColumn]) {
00602 marked[iColumn]=1;
00603 index[numberNonZero++]=iColumn;
00604 }
00605 array2[iColumn] -= value;
00606 }
00607 }
00608
00609 numberOriginal=numberNonZero;
00610 numberNonZero=0;
00611 for (i=0;i<numberOriginal;i++) {
00612 int iColumn = index[i];
00613 if (marked[iColumn]) {
00614 double value = array2[iColumn];
00615 array2[iColumn]=0.0;
00616 marked[iColumn]=0;
00617 if (fabs(value)>zeroTolerance) {
00618 array[numberNonZero]=value;
00619 index[numberNonZero++]=iColumn;
00620 }
00621 }
00622 }
00623 } else {
00624 for (i=0;i<numberOriginal;i++) {
00625 int iColumn = mark[i];
00626 index[i]=iColumn;
00627 array[iColumn]=markVector[iColumn];
00628 markVector[iColumn]=0.0;
00629 }
00630 numberNonZero=numberOriginal;
00631
00632 char * marked = (char *) markVector;
00633 for (i=0;i<numberOriginal;i++) {
00634 int iColumn = index[i];
00635 marked[iColumn]=0;
00636 }
00637 for (i=0;i<numberInRowArray;i++) {
00638 iRow = whichRow[i];
00639 double value = pi[iRow]*scalar;
00640 CoinBigIndex j;
00641 for (j=startPositive[iRow];j<startNegative[iRow];j++) {
00642 int iColumn = column[j];
00643 if (!marked[iColumn]) {
00644 marked[iColumn]=1;
00645 index[numberNonZero++]=iColumn;
00646 }
00647 array[iColumn] += value;
00648 }
00649 for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
00650 int iColumn = column[j];
00651 if (!marked[iColumn]) {
00652 marked[iColumn]=1;
00653 index[numberNonZero++]=iColumn;
00654 }
00655 array[iColumn] -= value;
00656 }
00657 }
00658
00659 numberOriginal=numberNonZero;
00660 numberNonZero=0;
00661 for (i=0;i<numberOriginal;i++) {
00662 int iColumn = index[i];
00663 marked[iColumn]=0;
00664 if (fabs(array[iColumn])>zeroTolerance) {
00665 index[numberNonZero++]=iColumn;
00666 } else {
00667 array[iColumn]=0.0;
00668 }
00669 }
00670 }
00671 } else if (numberInRowArray==2) {
00672
00673
00674 int iRow0 = whichRow[0];
00675 int iRow1 = whichRow[1];
00676 int j;
00677 if (packed) {
00678 double pi0 = pi[0];
00679 double pi1 = pi[1];
00680 if (startPositive[iRow0+1]-startPositive[iRow0]>
00681 startPositive[iRow1+1]-startPositive[iRow1]) {
00682 int temp = iRow0;
00683 iRow0 = iRow1;
00684 iRow1 = temp;
00685 pi0=pi1;
00686 pi1=pi[0];
00687 }
00688
00689 char * marked = (char *) (index+columnArray->capacity());
00690 int * lookup = y->getIndices();
00691 double value = pi0*scalar;
00692 for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
00693 int iColumn = column[j];
00694 array[numberNonZero] = value;
00695 marked[iColumn]=1;
00696 lookup[iColumn]=numberNonZero;
00697 index[numberNonZero++]=iColumn;
00698 }
00699 for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
00700 int iColumn = column[j];
00701 array[numberNonZero] = -value;
00702 marked[iColumn]=1;
00703 lookup[iColumn]=numberNonZero;
00704 index[numberNonZero++]=iColumn;
00705 }
00706 int numberOriginal = numberNonZero;
00707 value = pi1*scalar;
00708 for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
00709 int iColumn = column[j];
00710 if (marked[iColumn]) {
00711 int iLookup = lookup[iColumn];
00712 array[iLookup] += value;
00713 } else {
00714 if (fabs(value)>zeroTolerance) {
00715 array[numberNonZero] = value;
00716 index[numberNonZero++]=iColumn;
00717 }
00718 }
00719 }
00720 for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
00721 int iColumn = column[j];
00722 if (marked[iColumn]) {
00723 int iLookup = lookup[iColumn];
00724 array[iLookup] -= value;
00725 } else {
00726 if (fabs(value)>zeroTolerance) {
00727 array[numberNonZero] = -value;
00728 index[numberNonZero++]=iColumn;
00729 }
00730 }
00731 }
00732
00733 int nDelete=0;
00734 for (j=0;j<numberOriginal;j++) {
00735 int iColumn = index[j];
00736 marked[iColumn]=0;
00737 if (fabs(array[j])<=zeroTolerance)
00738 nDelete++;
00739 }
00740 if (nDelete) {
00741 numberOriginal=numberNonZero;
00742 numberNonZero=0;
00743 for (j=0;j<numberOriginal;j++) {
00744 int iColumn = index[j];
00745 double value = array[j];
00746 array[j]=0.0;
00747 if (fabs(value)>zeroTolerance) {
00748 array[numberNonZero]=value;
00749 index[numberNonZero++]=iColumn;
00750 }
00751 }
00752 }
00753 } else {
00754 if (startPositive[iRow0+1]-startPositive[iRow0]<
00755 startPositive[iRow1+1]-startPositive[iRow1]) {
00756 int temp = iRow0;
00757 iRow0 = iRow1;
00758 iRow1 = temp;
00759 }
00760 int numberOriginal;
00761 int i;
00762 numberNonZero=0;
00763 double value;
00764 value = pi[iRow0]*scalar;
00765 CoinBigIndex j;
00766 for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
00767 int iColumn = column[j];
00768 index[numberNonZero++]=iColumn;
00769 array[iColumn] = value;
00770 }
00771 for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
00772 int iColumn = column[j];
00773 index[numberNonZero++]=iColumn;
00774 array[iColumn] = -value;
00775 }
00776 value = pi[iRow1]*scalar;
00777 for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
00778 int iColumn = column[j];
00779 double value2= array[iColumn];
00780 if (value2) {
00781 value2 += value;
00782 } else {
00783 value2 = value;
00784 index[numberNonZero++]=iColumn;
00785 }
00786 array[iColumn] = value2;
00787 }
00788 for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
00789 int iColumn = column[j];
00790 double value2= array[iColumn];
00791 if (value2) {
00792 value2 -= value;
00793 } else {
00794 value2 = -value;
00795 index[numberNonZero++]=iColumn;
00796 }
00797 array[iColumn] = value2;
00798 }
00799
00800 numberOriginal=numberNonZero;
00801 numberNonZero=0;
00802 for (i=0;i<numberOriginal;i++) {
00803 int iColumn = index[i];
00804 if (fabs(array[iColumn])>zeroTolerance) {
00805 index[numberNonZero++]=iColumn;
00806 } else {
00807 array[iColumn]=0.0;
00808 }
00809 }
00810 }
00811 } else if (numberInRowArray==1) {
00812
00813 int iRow=rowArray->getIndices()[0];
00814 numberNonZero=0;
00815 double value;
00816 iRow = whichRow[0];
00817 CoinBigIndex j;
00818 if (packed) {
00819 value = pi[0]*scalar;
00820 if (fabs(value)>zeroTolerance) {
00821 for (j=startPositive[iRow];j<startNegative[iRow];j++) {
00822 int iColumn = column[j];
00823 array[numberNonZero] = value;
00824 index[numberNonZero++]=iColumn;
00825 }
00826 for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
00827 int iColumn = column[j];
00828 array[numberNonZero] = -value;
00829 index[numberNonZero++]=iColumn;
00830 }
00831 }
00832 } else {
00833 value = pi[iRow]*scalar;
00834 if (fabs(value)>zeroTolerance) {
00835 for (j=startPositive[iRow];j<startNegative[iRow];j++) {
00836 int iColumn = column[j];
00837 array[iColumn] = value;
00838 index[numberNonZero++]=iColumn;
00839 }
00840 for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
00841 int iColumn = column[j];
00842 array[iColumn] = -value;
00843 index[numberNonZero++]=iColumn;
00844 }
00845 }
00846 }
00847 }
00848 columnArray->setNumElements(numberNonZero);
00849 y->setNumElements(0);
00850 }
00851
00852
00853
00854 void
00855 ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex * model,
00856 const CoinIndexedVector * rowArray,
00857 const CoinIndexedVector * y,
00858 CoinIndexedVector * columnArray) const
00859 {
00860 columnArray->clear();
00861 double * pi = rowArray->denseVector();
00862 double * array = columnArray->denseVector();
00863 int jColumn;
00864 int numberToDo = y->getNumElements();
00865 const int * which = y->getIndices();
00866 assert (!rowArray->packedMode());
00867 columnArray->setPacked();
00868 for (jColumn=0;jColumn<numberToDo;jColumn++) {
00869 int iColumn = which[jColumn];
00870 double value = 0.0;
00871 CoinBigIndex j=startPositive_[iColumn];
00872 for (;j<startNegative_[iColumn];j++) {
00873 int iRow = indices_[j];
00874 value += pi[iRow];
00875 }
00876 for (;j<startPositive_[iColumn+1];j++) {
00877 int iRow = indices_[j];
00878 value -= pi[iRow];
00879 }
00880 array[jColumn]=value;
00881 }
00882 }
00883
00884
00885 CoinBigIndex
00886 ClpPlusMinusOneMatrix::numberInBasis(const int * columnIsBasic) const
00887 {
00888 int i;
00889 CoinBigIndex numberElements=0;
00890 assert (columnOrdered_);
00891 for (i=0;i<numberColumns_;i++) {
00892 if (columnIsBasic[i]>=0)
00893 numberElements += startPositive_[i+1]-startPositive_[i];
00894 }
00895 return numberElements;
00896 }
00897
00898 CoinBigIndex
00899 ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
00900 const int * columnIsBasic, int & numberBasic,
00901 int * indexRowU, int * indexColumnU,
00902 double * elementU) const
00903 {
00904 #ifdef CLPDEBUG
00905 const double * rowScale = model->rowScale();
00906 assert (!rowScale);
00907 #endif
00908 int i;
00909 CoinBigIndex numberElements=0;
00910 assert (columnOrdered_);
00911 for (i=0;i<numberColumns_;i++) {
00912 if (columnIsBasic[i]>=0) {
00913 CoinBigIndex j=startPositive_[i];
00914 for (;j<startNegative_[i];j++) {
00915 int iRow = indices_[j];
00916 indexRowU[numberElements]=iRow;
00917 indexColumnU[numberElements]=numberBasic;
00918 elementU[numberElements++]=1.0;
00919 }
00920 for (;j<startPositive_[i+1];j++) {
00921 int iRow = indices_[j];
00922 indexRowU[numberElements]=iRow;
00923 indexColumnU[numberElements]=numberBasic;
00924 elementU[numberElements++]=-1.0;
00925 }
00926 numberBasic++;
00927 }
00928 }
00929 return numberElements;
00930 }
00931
00932
00933 CoinBigIndex
00934 ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
00935 const int * whichColumn,
00936 int numberBasic,
00937 int numberColumnBasic,
00938 int * indexRowU, int * indexColumnU,
00939 double * elementU) const
00940 {
00941 int i;
00942 CoinBigIndex numberElements=0;
00943 if (elementU!=NULL) {
00944 assert (columnOrdered_);
00945 for (i=0;i<numberColumnBasic;i++) {
00946 int iColumn = whichColumn[i];
00947 CoinBigIndex j=startPositive_[iColumn];
00948 for (;j<startNegative_[iColumn];j++) {
00949 int iRow = indices_[j];
00950 indexRowU[numberElements]=iRow;
00951 indexColumnU[numberElements]=numberBasic;
00952 elementU[numberElements++]=1.0;
00953 }
00954 for (;j<startPositive_[iColumn+1];j++) {
00955 int iRow = indices_[j];
00956 indexRowU[numberElements]=iRow;
00957 indexColumnU[numberElements]=numberBasic;
00958 elementU[numberElements++]=-1.0;
00959 }
00960 numberBasic++;
00961 }
00962 } else {
00963 for (i=0;i<numberColumnBasic;i++) {
00964 int iColumn = whichColumn[i];
00965 numberElements += startPositive_[iColumn+1]-startPositive_[iColumn];
00966 }
00967 }
00968 return numberElements;
00969 }
00970
00971
00972 void
00973 ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model,
00974 CoinIndexedVector * rowArray,
00975 int iColumn) const
00976 {
00977 CoinBigIndex j=startPositive_[iColumn];
00978 for (;j<startNegative_[iColumn];j++) {
00979 int iRow = indices_[j];
00980 rowArray->add(iRow,1.0);
00981 }
00982 for (;j<startPositive_[iColumn+1];j++) {
00983 int iRow = indices_[j];
00984 rowArray->add(iRow,-1.0);
00985 }
00986 }
00987
00988
00989
00990
00991 void
00992 ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * model,
00993 CoinIndexedVector * rowArray,
00994 int iColumn) const
00995 {
00996 int * index = rowArray->getIndices();
00997 double * array = rowArray->denseVector();
00998 int number = 0;
00999 CoinBigIndex j=startPositive_[iColumn];
01000 for (;j<startNegative_[iColumn];j++) {
01001 int iRow = indices_[j];
01002 array[number]=1.0;
01003 index[number++]=iRow;
01004 }
01005 for (;j<startPositive_[iColumn+1];j++) {
01006 int iRow = indices_[j];
01007 array[number]=-1.0;
01008 index[number++]=iRow;
01009 }
01010 rowArray->setNumElements(number);
01011 rowArray->setPackedMode(true);
01012 }
01013
01014
01015 void
01016 ClpPlusMinusOneMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
01017 int iColumn, double multiplier) const
01018 {
01019 CoinBigIndex j=startPositive_[iColumn];
01020 for (;j<startNegative_[iColumn];j++) {
01021 int iRow = indices_[j];
01022 rowArray->quickAdd(iRow,multiplier);
01023 }
01024 for (;j<startPositive_[iColumn+1];j++) {
01025 int iRow = indices_[j];
01026 rowArray->quickAdd(iRow,-multiplier);
01027 }
01028 }
01029
01030
01031 CoinPackedMatrix *
01032 ClpPlusMinusOneMatrix::getPackedMatrix() const
01033 {
01034 int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
01035 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
01036 return new CoinPackedMatrix(columnOrdered_,numberMinor,numberMajor,
01037 getNumElements(),
01038 getElements(),indices_,
01039 startPositive_,getVectorLengths());
01040
01041 }
01042
01043
01044
01045
01046 const double *
01047 ClpPlusMinusOneMatrix::getElements() const
01048 {
01049 if (!elements_) {
01050 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
01051 int numberElements = startPositive_[numberMajor];
01052 elements_ = new double [numberElements];
01053 CoinBigIndex j=0;
01054 int i;
01055 for (i=0;i<numberMajor;i++) {
01056 for (;j<startNegative_[i];j++) {
01057 elements_[j]=1.0;
01058 }
01059 for (;j<startPositive_[i+1];j++) {
01060 elements_[j]=-1.0;
01061 }
01062 }
01063 }
01064 return elements_;
01065 }
01066
01067 const CoinBigIndex *
01068 ClpPlusMinusOneMatrix::getVectorStarts() const
01069 {
01070 return startPositive_;
01071 }
01072
01073 const int *
01074 ClpPlusMinusOneMatrix::getVectorLengths() const
01075 {
01076 if (!lengths_) {
01077 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
01078 lengths_ = new int [numberMajor];
01079 int i;
01080 for (i=0;i<numberMajor;i++) {
01081 lengths_[i]=startPositive_[i+1]-startPositive_[i];
01082 }
01083 }
01084 return lengths_;
01085 }
01086
01087 void
01088 ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel)
01089 {
01090 int iColumn;
01091 int newSize=startPositive_[numberColumns_];;
01092 int numberBad=0;
01093
01094 int * which = new int[numberColumns_];
01095 memset(which,0,numberColumns_*sizeof(int));
01096 int nDuplicate=0;
01097 for (iColumn=0;iColumn<numDel;iColumn++) {
01098 int jColumn = indDel[iColumn];
01099 if (jColumn<0||jColumn>=numberColumns_) {
01100 numberBad++;
01101 } else {
01102 newSize -= startPositive_[jColumn+1]-startPositive_[jColumn];
01103 if (which[jColumn])
01104 nDuplicate++;
01105 else
01106 which[jColumn]=1;
01107 }
01108 }
01109 if (numberBad)
01110 throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
01111 int newNumber = numberColumns_-numDel+nDuplicate;
01112
01113 delete [] lengths_;
01114 lengths_=NULL;
01115 delete [] elements_;
01116 elements_= NULL;
01117 int * newPositive = new int [newNumber+1];
01118 int * newNegative = new int [newNumber];
01119 int * newIndices = new int [newSize];
01120 newNumber=0;
01121 newSize=0;
01122 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
01123 if (!which[iColumn]) {
01124 int start,end;
01125 int i;
01126 start = startPositive_[iColumn];
01127 end=startNegative_[iColumn];
01128 newPositive[newNumber]=newSize;
01129 for (i=start;i<end;i++)
01130 newIndices[newSize++]=indices_[i];
01131 start = startNegative_[iColumn];
01132 end=startPositive_[iColumn+1];
01133 newNegative[newNumber++]=newSize;
01134 for (i=start;i<end;i++)
01135 newIndices[newSize++]=indices_[i];
01136 }
01137 }
01138 newPositive[newNumber]=newSize;
01139 delete [] which;
01140 delete [] startPositive_;
01141 startPositive_= newPositive;
01142 delete [] startNegative_;
01143 startNegative_= newNegative;
01144 delete [] indices_;
01145 indices_= newIndices;
01146 numberColumns_ = newNumber;
01147 }
01148
01149 void
01150 ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel)
01151 {
01152 int iRow;
01153 int numberBad=0;
01154
01155 int * which = new int[numberRows_];
01156 memset(which,0,numberRows_*sizeof(int));
01157 int nDuplicate=0;
01158 for (iRow=0;iRow<numDel;iRow++) {
01159 int jRow = indDel[iRow];
01160 if (jRow<0||jRow>=numberRows_) {
01161 numberBad++;
01162 } else {
01163 if (which[jRow])
01164 nDuplicate++;
01165 else
01166 which[jRow]=1;
01167 }
01168 }
01169 if (numberBad)
01170 throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
01171 int iElement;
01172 int numberElements=startPositive_[numberColumns_];
01173 int newSize=0;
01174 for (iElement=0;iElement<numberElements;iElement++) {
01175 iRow = indices_[iElement];
01176 if (!which[iRow])
01177 newSize++;
01178 }
01179 int newNumber = numberRows_-numDel+nDuplicate;
01180
01181 delete [] lengths_;
01182 lengths_=NULL;
01183 delete [] elements_;
01184 elements_= NULL;
01185 int * newIndices = new int [newSize];
01186 newSize=0;
01187 int iColumn;
01188 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
01189 int start,end;
01190 int i;
01191 start = startPositive_[iColumn];
01192 end=startNegative_[iColumn];
01193 startPositive_[newNumber]=newSize;
01194 for (i=start;i<end;i++) {
01195 iRow = indices_[i];
01196 if (!which[iRow])
01197 newIndices[newSize++]=iRow;
01198 }
01199 start = startNegative_[iColumn];
01200 end=startPositive_[iColumn+1];
01201 startNegative_[newNumber]=newSize;
01202 for (i=start;i<end;i++) {
01203 iRow = indices_[i];
01204 if (!which[iRow])
01205 newIndices[newSize++]=iRow;
01206 }
01207 }
01208 startPositive_[numberColumns_]=newSize;
01209 delete [] which;
01210 delete [] indices_;
01211 indices_= newIndices;
01212 numberRows_ = newNumber;
01213 }
01214 bool
01215 ClpPlusMinusOneMatrix::isColOrdered() const
01216 {
01217 return columnOrdered_;
01218 }
01219
01220 CoinBigIndex
01221 ClpPlusMinusOneMatrix::getNumElements() const
01222 {
01223 int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
01224 if (startPositive_)
01225 return startPositive_[numberMajor];
01226 else
01227 return 0;
01228 }
01229
01230 void
01231 ClpPlusMinusOneMatrix::passInCopy(int numberRows, int numberColumns,
01232 bool columnOrdered, int * indices,
01233 int * startPositive, int * startNegative)
01234 {
01235 columnOrdered_=columnOrdered;
01236 startPositive_ = startPositive;
01237 startNegative_ = startNegative;
01238 indices_ = indices;
01239 numberRows_=numberRows;
01240 numberColumns_=numberColumns;
01241 }
01242
01243
01244
01245
01246 CoinBigIndex *
01247 ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
01248 {
01249 int numberRows = model->numberRows();
01250 int numberColumns =model->numberColumns();
01251 int number = numberRows+numberColumns;
01252 CoinBigIndex * weights = new CoinBigIndex[number];
01253 int i;
01254 for (i=0;i<numberColumns;i++) {
01255 CoinBigIndex j;
01256 CoinBigIndex count=0;
01257 for (j=startPositive_[i];j<startPositive_[i+1];j++) {
01258 int iRow=indices_[j];
01259 count += inputWeights[iRow];
01260 }
01261 weights[i]=count;
01262 }
01263 for (i=0;i<numberRows;i++) {
01264 weights[i+numberColumns]=inputWeights[i];
01265 }
01266 return weights;
01267 }
01268
01269 void
01270 ClpPlusMinusOneMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
01271 {
01272 int iColumn;
01273 int size=0;
01274 int numberBad=0;
01275 for (iColumn=0;iColumn<number;iColumn++) {
01276 int n=columns[iColumn]->getNumElements();
01277 const double * element = columns[iColumn]->getElements();
01278 size += n;
01279 int i;
01280 for (i=0;i<n;i++) {
01281 if (fabs(element[i])!=1.0)
01282 numberBad++;
01283 }
01284 }
01285 if (numberBad)
01286 throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix");
01287
01288 delete [] lengths_;
01289 lengths_=NULL;
01290 delete [] elements_;
01291 elements_= NULL;
01292 int numberNow = startPositive_[numberColumns_];
01293 int * temp;
01294 temp = new int [numberColumns_+1+number];
01295 memcpy(temp,startPositive_,(numberColumns_+1)*sizeof(int));
01296 delete [] startPositive_;
01297 startPositive_= temp;
01298 temp = new int [numberColumns_+number];
01299 memcpy(temp,startNegative_,numberColumns_*sizeof(int));
01300 delete [] startNegative_;
01301 startNegative_= temp;
01302 temp = new int [numberNow+size];
01303 memcpy(temp,indices_,numberNow*sizeof(int));
01304 delete [] indices_;
01305 indices_= temp;
01306
01307 size=numberNow;
01308 for (iColumn=0;iColumn<number;iColumn++) {
01309 int n=columns[iColumn]->getNumElements();
01310 const int * row = columns[iColumn]->getIndices();
01311 const double * element = columns[iColumn]->getElements();
01312 int i;
01313 for (i=0;i<n;i++) {
01314 if (element[i]==1.0)
01315 indices_[size++]=row[i];
01316 }
01317 startNegative_[iColumn+numberColumns_]=size;
01318 for (i=0;i<n;i++) {
01319 if (element[i]==-1.0)
01320 indices_[size++]=row[i];
01321 }
01322 startPositive_[iColumn+numberColumns_+1]=size;
01323 }
01324
01325 numberColumns_ += number;
01326 }
01327
01328 void
01329 ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
01330 {
01331
01332 int * countPositive = new int [numberColumns_+1];
01333 memset(countPositive,0,numberColumns_*sizeof(int));
01334 int * countNegative = new int [numberColumns_];
01335 memset(countNegative,0,numberColumns_*sizeof(int));
01336 int iRow;
01337 int size=0;
01338 int numberBad=0;
01339 for (iRow=0;iRow<number;iRow++) {
01340 int n=rows[iRow]->getNumElements();
01341 const int * column = rows[iRow]->getIndices();
01342 const double * element = rows[iRow]->getElements();
01343 size += n;
01344 int i;
01345 for (i=0;i<n;i++) {
01346 int iColumn = column[i];
01347 if (element[i]==1.0)
01348 countPositive[iColumn]++;
01349 else if (element[i]==-1.0)
01350 countNegative[iColumn]++;
01351 else
01352 numberBad++;
01353 }
01354 }
01355 if (numberBad)
01356 throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix");
01357
01358 delete [] lengths_;
01359 lengths_=NULL;
01360 delete [] elements_;
01361 elements_= NULL;
01362 int numberNow = startPositive_[numberColumns_];
01363 int * newIndices = new int [numberNow+size];
01364
01365
01366 int iColumn;
01367 int numberAdded=0;
01368 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
01369 int n,now,move;
01370 now = startPositive_[iColumn];
01371 move = startNegative_[iColumn]-now;
01372 n = countPositive[iColumn];
01373 startPositive_[iColumn] += numberAdded;
01374 memcpy(indices_+now,newIndices+startPositive_[iColumn],move*sizeof(int));
01375 countPositive[iColumn]= startNegative_[iColumn]+numberAdded;
01376 numberAdded += n;
01377 now = startNegative_[iColumn];
01378 move = startPositive_[iColumn+1]-now;
01379 n = countNegative[iColumn];
01380 startNegative_[iColumn] += numberAdded;
01381 memcpy(indices_+now,newIndices+startNegative_[iColumn],move*sizeof(int));
01382 countNegative[iColumn]= startPositive_[iColumn+1]+numberAdded;
01383 numberAdded += n;
01384 }
01385 delete [] indices_;
01386 indices_ = newIndices;
01387 startPositive_[numberColumns_] += numberAdded;
01388
01389 for (iRow=0;iRow<number;iRow++) {
01390 int newRow = numberRows_+iRow;
01391 int n=rows[iRow]->getNumElements();
01392 const int * column = rows[iRow]->getIndices();
01393 const double * element = rows[iRow]->getElements();
01394 int i;
01395 for (i=0;i<n;i++) {
01396 int iColumn = column[i];
01397 int put;
01398 if (element[i]==1.0) {
01399 put = countPositive[iColumn];
01400 countPositive[iColumn] = put+1;
01401 } else {
01402 put = countNegative[iColumn];
01403 countNegative[iColumn] = put+1;
01404 }
01405 indices_[put]=newRow;
01406 }
01407 }
01408 delete [] countPositive;
01409 delete [] countNegative;
01410 numberRows_ += number;
01411 }
01412
01413
01414
01415 void
01416 ClpPlusMinusOneMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
01417 double & smallestPositive, double & largestPositive)
01418 {
01419 int iColumn;
01420 bool plusOne=false;
01421 bool minusOne=false;
01422 for (iColumn=0;iColumn<numberColumns_;iColumn++) {
01423 if (startNegative_[iColumn]>startPositive_[iColumn])
01424 plusOne=true;
01425 if (startPositive_[iColumn+1]>startNegative_[iColumn])
01426 minusOne=true;
01427 }
01428 if (minusOne) {
01429 smallestNegative=-1.0;
01430 largestNegative=-1.0;
01431 } else {
01432 smallestNegative=0.0;
01433 largestNegative=0.0;
01434 }
01435 if (plusOne) {
01436 smallestPositive=1.0;
01437 largestPositive=1.0;
01438 } else {
01439 smallestPositive=0.0;
01440 largestPositive=0.0;
01441 }
01442 }
01443
01444 bool
01445 ClpPlusMinusOneMatrix::canDoPartialPricing() const
01446 {
01447 return true;
01448 }
01449
01450 void
01451 ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, int start, int end,
01452 int & bestSequence, int & numberWanted)
01453 {
01454 int j;
01455 double tolerance=model->currentDualTolerance();
01456 double * reducedCost = model->djRegion();
01457 const double * duals = model->dualRowSolution();
01458 const double * cost = model->costRegion();
01459 double bestDj;
01460 if (bestSequence>=0)
01461 bestDj = fabs(reducedCost[bestSequence]);
01462 else
01463 bestDj=tolerance;
01464 int sequenceOut = model->sequenceOut();
01465 int saveSequence = bestSequence;
01466 int iSequence;
01467 for (iSequence=start;iSequence<end;iSequence++) {
01468 if (iSequence!=sequenceOut) {
01469 double value;
01470 ClpSimplex::Status status = model->getStatus(iSequence);
01471
01472 switch(status) {
01473
01474 case ClpSimplex::basic:
01475 case ClpSimplex::isFixed:
01476 break;
01477 case ClpSimplex::isFree:
01478 case ClpSimplex::superBasic:
01479 value=cost[iSequence];
01480 j=startPositive_[iSequence];
01481 for (;j<startNegative_[iSequence];j++) {
01482 int iRow = indices_[j];
01483 value -= duals[iRow];
01484 }
01485 for (;j<startPositive_[iSequence+1];j++) {
01486 int iRow = indices_[j];
01487 value += duals[iRow];
01488 }
01489 value = fabs(value);
01490 if (value>FREE_ACCEPT*tolerance) {
01491 numberWanted--;
01492
01493 value *= FREE_BIAS;
01494 if (value>bestDj) {
01495
01496 if (!model->flagged(iSequence)) {
01497 bestDj=value;
01498 bestSequence = iSequence;
01499 } else {
01500
01501 numberWanted++;
01502 }
01503 }
01504 }
01505 break;
01506 case ClpSimplex::atUpperBound:
01507 value=cost[iSequence];
01508 j=startPositive_[iSequence];
01509 for (;j<startNegative_[iSequence];j++) {
01510 int iRow = indices_[j];
01511 value -= duals[iRow];
01512 }
01513 for (;j<startPositive_[iSequence+1];j++) {
01514 int iRow = indices_[j];
01515 value += duals[iRow];
01516 }
01517 if (value>tolerance) {
01518 numberWanted--;
01519 if (value>bestDj) {
01520
01521 if (!model->flagged(iSequence)) {
01522 bestDj=value;
01523 bestSequence = iSequence;
01524 } else {
01525
01526 numberWanted++;
01527 }
01528 }
01529 }
01530 break;
01531 case ClpSimplex::atLowerBound:
01532 value=cost[iSequence];
01533 j=startPositive_[iSequence];
01534 for (;j<startNegative_[iSequence];j++) {
01535 int iRow = indices_[j];
01536 value -= duals[iRow];
01537 }
01538 for (;j<startPositive_[iSequence+1];j++) {
01539 int iRow = indices_[j];
01540 value += duals[iRow];
01541 }
01542 value = -value;
01543 if (value>tolerance) {
01544 numberWanted--;
01545 if (value>bestDj) {
01546
01547 if (!model->flagged(iSequence)) {
01548 bestDj=value;
01549 bestSequence = iSequence;
01550 } else {
01551
01552 numberWanted++;
01553 }
01554 }
01555 }
01556 break;
01557 }
01558 }
01559 if (!numberWanted)
01560 break;
01561 }
01562 if (bestSequence!=saveSequence) {
01563
01564 double value=cost[bestSequence];
01565 j=startPositive_[bestSequence];
01566 for (;j<startNegative_[bestSequence];j++) {
01567 int iRow = indices_[j];
01568 value -= duals[iRow];
01569 }
01570 for (;j<startPositive_[bestSequence+1];j++) {
01571 int iRow = indices_[j];
01572 value += duals[iRow];
01573 }
01574 reducedCost[bestSequence] = value;
01575 }
01576 }