Changeset 2385 for trunk/Clp/src/ClpNetworkMatrix.cpp
 Timestamp:
 Jan 6, 2019 2:43:06 PM (3 months ago)
 File:

 1 edited
Legend:
 Unmodified
 Added
 Removed

trunk/Clp/src/ClpNetworkMatrix.cpp
r2271 r2385 3 3 // Corporation and others. All Rights Reserved. 4 4 // This code is licensed under the terms of the Eclipse Public License (EPL). 5 6 5 7 6 #include <cstdio> … … 28 27 // Default Constructor 29 28 // 30 ClpNetworkMatrix::ClpNetworkMatrix 31 32 { 33 34 35 36 37 38 39 29 ClpNetworkMatrix::ClpNetworkMatrix() 30 : ClpMatrixBase() 31 { 32 setType(11); 33 matrix_ = NULL; 34 lengths_ = NULL; 35 indices_ = NULL; 36 numberRows_ = 0; 37 numberColumns_ = 0; 38 trueNetwork_ = false; 40 39 } 41 40 42 41 /* Constructor from two arrays */ 43 ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int * head, 44 const int * tail) 45 : ClpMatrixBase() 46 { 47 setType(11); 48 matrix_ = NULL; 49 lengths_ = NULL; 50 indices_ = new int[2*numberColumns];; 51 numberRows_ = 1; 52 numberColumns_ = numberColumns; 53 trueNetwork_ = true; 54 int iColumn; 55 CoinBigIndex j = 0; 56 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 57 int iRow = head[iColumn]; 58 numberRows_ = CoinMax(numberRows_, iRow); 59 indices_[j] = iRow; 60 iRow = tail[iColumn]; 61 numberRows_ = CoinMax(numberRows_, iRow); 62 indices_[j+1] = iRow; 63 } 64 numberRows_++; 42 ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int *head, 43 const int *tail) 44 : ClpMatrixBase() 45 { 46 setType(11); 47 matrix_ = NULL; 48 lengths_ = NULL; 49 indices_ = new int[2 * numberColumns]; 50 ; 51 numberRows_ = 1; 52 numberColumns_ = numberColumns; 53 trueNetwork_ = true; 54 int iColumn; 55 CoinBigIndex j = 0; 56 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 57 int iRow = head[iColumn]; 58 numberRows_ = CoinMax(numberRows_, iRow); 59 indices_[j] = iRow; 60 iRow = tail[iColumn]; 61 numberRows_ = CoinMax(numberRows_, iRow); 62 indices_[j + 1] = iRow; 63 } 64 numberRows_++; 65 65 } 66 66 // 67 67 // Copy constructor 68 68 // 69 ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix &rhs)70 71 { 72 73 74 75 76 77 78 79 indices_ = new int [ 2*numberColumns_];80 81 82 83 84 85 86 87 88 } 89 90 ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix &rhs)91 92 { 93 94 95 96 97 98 assert(rhs.isColOrdered());99 100 const int *row = rhs.getIndices();101 const CoinBigIndex *columnStart = rhs.getVectorStarts();102 const int *columnLength = rhs.getVectorLengths();103 const double *elementByColumn = rhs.getElements();104 105 106 107 indices_ = new int[2*numberColumns_];108 109 110 111 112 113 114 115 116 indices_[j+1] = 1;117 118 119 120 121 122 123 124 125 indices_[j+1] = iRow;126 127 indices_[j+1] = 1;128 129 130 131 132 133 134 135 136 137 138 if (fabs(elementByColumn[k+1] + 1.0) < 1.0e10) {139 140 141 indices_[j+1] = iRow;142 iRow = row[k+1];143 144 145 146 147 148 149 if (fabs(elementByColumn[k+1]  1.0) < 1.0e10) {150 151 152 153 iRow = row[k+1];154 155 indices_[j+1] = iRow;156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 delete[] indices_;173 174 175 176 177 178 179 numberRows_++; // correct180 181 69 ClpNetworkMatrix::ClpNetworkMatrix(const ClpNetworkMatrix &rhs) 70 : ClpMatrixBase(rhs) 71 { 72 matrix_ = NULL; 73 lengths_ = NULL; 74 indices_ = NULL; 75 numberRows_ = rhs.numberRows_; 76 numberColumns_ = rhs.numberColumns_; 77 trueNetwork_ = rhs.trueNetwork_; 78 if (numberColumns_) { 79 indices_ = new int[2 * numberColumns_]; 80 CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_); 81 } 82 int numberRows = getNumRows(); 83 if (rhs.rhsOffset_ && numberRows) { 84 rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows); 85 } else { 86 rhsOffset_ = NULL; 87 } 88 } 89 90 ClpNetworkMatrix::ClpNetworkMatrix(const CoinPackedMatrix &rhs) 91 : ClpMatrixBase() 92 { 93 setType(11); 94 matrix_ = NULL; 95 lengths_ = NULL; 96 indices_ = NULL; 97 int iColumn; 98 assert(rhs.isColOrdered()); 99 // get matrix data pointers 100 const int *row = rhs.getIndices(); 101 const CoinBigIndex *columnStart = rhs.getVectorStarts(); 102 const int *columnLength = rhs.getVectorLengths(); 103 const double *elementByColumn = rhs.getElements(); 104 numberColumns_ = rhs.getNumCols(); 105 int goodNetwork = 1; 106 numberRows_ = 1; 107 indices_ = new int[2 * numberColumns_]; 108 CoinBigIndex j = 0; 109 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 110 CoinBigIndex k = columnStart[iColumn]; 111 int iRow; 112 switch (columnLength[iColumn]) { 113 case 0: 114 goodNetwork = 1; // not classic network 115 indices_[j] = 1; 116 indices_[j + 1] = 1; 117 break; 118 119 case 1: 120 goodNetwork = 1; // not classic network 121 if (fabs(elementByColumn[k]  1.0) < 1.0e10) { 122 indices_[j] = 1; 123 iRow = row[k]; 124 numberRows_ = CoinMax(numberRows_, iRow); 125 indices_[j + 1] = iRow; 126 } else if (fabs(elementByColumn[k] + 1.0) < 1.0e10) { 127 indices_[j + 1] = 1; 128 iRow = row[k]; 129 numberRows_ = CoinMax(numberRows_, iRow); 130 indices_[j] = iRow; 131 } else { 132 goodNetwork = 0; // not a network 133 } 134 break; 135 136 case 2: 137 if (fabs(elementByColumn[k]  1.0) < 1.0e10) { 138 if (fabs(elementByColumn[k + 1] + 1.0) < 1.0e10) { 139 iRow = row[k]; 140 numberRows_ = CoinMax(numberRows_, iRow); 141 indices_[j + 1] = iRow; 142 iRow = row[k + 1]; 143 numberRows_ = CoinMax(numberRows_, iRow); 144 indices_[j] = iRow; 145 } else { 146 goodNetwork = 0; // not a network 147 } 148 } else if (fabs(elementByColumn[k] + 1.0) < 1.0e10) { 149 if (fabs(elementByColumn[k + 1]  1.0) < 1.0e10) { 150 iRow = row[k]; 151 numberRows_ = CoinMax(numberRows_, iRow); 152 indices_[j] = iRow; 153 iRow = row[k + 1]; 154 numberRows_ = CoinMax(numberRows_, iRow); 155 indices_[j + 1] = iRow; 156 } else { 157 goodNetwork = 0; // not a network 158 } 159 } else { 160 goodNetwork = 0; // not a network 161 } 162 break; 163 164 default: 165 goodNetwork = 0; // not a network 166 break; 167 } 168 if (!goodNetwork) 169 break; 170 } 171 if (!goodNetwork) { 172 delete[] indices_; 173 // put in message 174 printf("Not a network  can test if indices_ null\n"); 175 indices_ = NULL; 176 numberRows_ = 0; 177 numberColumns_ = 0; 178 } else { 179 numberRows_++; // correct 180 trueNetwork_ = goodNetwork > 0; 181 } 182 182 } 183 183 … … 185 185 // Destructor 186 186 // 187 ClpNetworkMatrix::~ClpNetworkMatrix 188 { 189 190 delete[] lengths_;191 delete[] indices_;187 ClpNetworkMatrix::~ClpNetworkMatrix() 188 { 189 delete matrix_; 190 delete[] lengths_; 191 delete[] indices_; 192 192 } 193 193 … … 196 196 // 197 197 ClpNetworkMatrix & 198 ClpNetworkMatrix::operator=(const ClpNetworkMatrix &rhs)199 { 200 201 202 203 delete[] lengths_;204 delete[] indices_;205 206 207 208 209 210 211 212 indices_ = new int [ 2*numberColumns_];213 214 215 216 198 ClpNetworkMatrix::operator=(const ClpNetworkMatrix &rhs) 199 { 200 if (this != &rhs) { 201 ClpMatrixBase::operator=(rhs); 202 delete matrix_; 203 delete[] lengths_; 204 delete[] indices_; 205 matrix_ = NULL; 206 lengths_ = NULL; 207 indices_ = NULL; 208 numberRows_ = rhs.numberRows_; 209 numberColumns_ = rhs.numberColumns_; 210 trueNetwork_ = rhs.trueNetwork_; 211 if (numberColumns_) { 212 indices_ = new int[2 * numberColumns_]; 213 CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_); 214 } 215 } 216 return *this; 217 217 } 218 218 // 219 219 // Clone 220 220 // 221 ClpMatrixBase * 222 { 223 221 ClpMatrixBase *ClpNetworkMatrix::clone() const 222 { 223 return new ClpNetworkMatrix(*this); 224 224 } 225 225 … … 228 228 ClpNetworkMatrix::reverseOrderedCopy() const 229 229 { 230 231 CoinBigIndex * tempP = new CoinBigIndex[numberRows_];232 CoinBigIndex * tempN = new CoinBigIndex[numberRows_];233 234 235 236 237 238 239 240 iRow = indices_[j+1];241 242 243 int * newIndices = new int [2*numberColumns_];244 CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];245 CoinBigIndex *newN = new CoinBigIndex[numberRows_];246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 iRow = indices_[j+1];265 266 267 268 269 delete[] tempP;270 delete[] tempN;271 ClpPlusMinusOneMatrix *newCopy = new ClpPlusMinusOneMatrix();272 273 false,newIndices, newP, newN);274 230 // count number in each row 231 CoinBigIndex *tempP = new CoinBigIndex[numberRows_]; 232 CoinBigIndex *tempN = new CoinBigIndex[numberRows_]; 233 memset(tempP, 0, numberRows_ * sizeof(CoinBigIndex)); 234 memset(tempN, 0, numberRows_ * sizeof(CoinBigIndex)); 235 CoinBigIndex j = 0; 236 int i; 237 for (i = 0; i < numberColumns_; i++, j += 2) { 238 int iRow = indices_[j]; 239 tempN[iRow]++; 240 iRow = indices_[j + 1]; 241 tempP[iRow]++; 242 } 243 int *newIndices = new int[2 * numberColumns_]; 244 CoinBigIndex *newP = new CoinBigIndex[numberRows_ + 1]; 245 CoinBigIndex *newN = new CoinBigIndex[numberRows_]; 246 int iRow; 247 j = 0; 248 // do starts 249 for (iRow = 0; iRow < numberRows_; iRow++) { 250 newP[iRow] = j; 251 j += tempP[iRow]; 252 tempP[iRow] = newP[iRow]; 253 newN[iRow] = j; 254 j += tempN[iRow]; 255 tempN[iRow] = newN[iRow]; 256 } 257 newP[numberRows_] = j; 258 j = 0; 259 for (i = 0; i < numberColumns_; i++, j += 2) { 260 int iRow = indices_[j]; 261 CoinBigIndex put = tempN[iRow]; 262 newIndices[put++] = i; 263 tempN[iRow] = put; 264 iRow = indices_[j + 1]; 265 put = tempP[iRow]; 266 newIndices[put++] = i; 267 tempP[iRow] = put; 268 } 269 delete[] tempP; 270 delete[] tempN; 271 ClpPlusMinusOneMatrix *newCopy = new ClpPlusMinusOneMatrix(); 272 newCopy>passInCopy(numberRows_, numberColumns_, 273 false, newIndices, newP, newN); 274 return newCopy; 275 275 } 276 276 //unscaled versions 277 void 278 ClpNetworkMatrix::times(double scalar, 279 const double * x, double * y) const 280 { 281 int iColumn; 282 CoinBigIndex j = 0; 283 if (trueNetwork_) { 284 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 285 double value = scalar * x[iColumn]; 286 if (value) { 287 int iRowM = indices_[j]; 288 int iRowP = indices_[j+1]; 289 y[iRowM] = value; 290 y[iRowP] += value; 291 } 292 } 293 } else { 294 // skip negative rows 295 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 296 double value = scalar * x[iColumn]; 297 if (value) { 298 int iRowM = indices_[j]; 299 int iRowP = indices_[j+1]; 300 if (iRowM >= 0) 301 y[iRowM] = value; 302 if (iRowP >= 0) 303 y[iRowP] += value; 304 } 305 } 306 } 307 } 308 void 309 ClpNetworkMatrix::transposeTimes(double scalar, 310 const double * x, double * y) const 311 { 312 int iColumn; 313 CoinBigIndex j = 0; 314 if (trueNetwork_) { 315 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 316 double value = y[iColumn]; 317 int iRowM = indices_[j]; 318 int iRowP = indices_[j+1]; 319 value = scalar * x[iRowM]; 320 value += scalar * x[iRowP]; 321 y[iColumn] = value; 322 } 323 } else { 324 // skip negative rows 325 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 326 double value = y[iColumn]; 327 int iRowM = indices_[j]; 328 int iRowP = indices_[j+1]; 329 if (iRowM >= 0) 330 value = scalar * x[iRowM]; 331 if (iRowP >= 0) 332 value += scalar * x[iRowP]; 333 y[iColumn] = value; 334 } 335 } 336 } 337 void 338 ClpNetworkMatrix::times(double scalar, 339 const double * x, double * y, 340 const double * /*rowScale*/, 341 const double * /*columnScale*/) const 342 { 343 // we know it is not scaled 344 times(scalar, x, y); 345 } 346 void 347 ClpNetworkMatrix::transposeTimes( double scalar, 348 const double * x, double * y, 349 const double * /*rowScale*/, 350 const double * /*columnScale*/, 351 double * /*spare*/) const 352 { 353 // we know it is not scaled 354 transposeTimes(scalar, x, y); 277 void ClpNetworkMatrix::times(double scalar, 278 const double *x, double *y) const 279 { 280 int iColumn; 281 CoinBigIndex j = 0; 282 if (trueNetwork_) { 283 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 284 double value = scalar * x[iColumn]; 285 if (value) { 286 int iRowM = indices_[j]; 287 int iRowP = indices_[j + 1]; 288 y[iRowM] = value; 289 y[iRowP] += value; 290 } 291 } 292 } else { 293 // skip negative rows 294 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 295 double value = scalar * x[iColumn]; 296 if (value) { 297 int iRowM = indices_[j]; 298 int iRowP = indices_[j + 1]; 299 if (iRowM >= 0) 300 y[iRowM] = value; 301 if (iRowP >= 0) 302 y[iRowP] += value; 303 } 304 } 305 } 306 } 307 void ClpNetworkMatrix::transposeTimes(double scalar, 308 const double *x, double *y) const 309 { 310 int iColumn; 311 CoinBigIndex j = 0; 312 if (trueNetwork_) { 313 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 314 double value = y[iColumn]; 315 int iRowM = indices_[j]; 316 int iRowP = indices_[j + 1]; 317 value = scalar * x[iRowM]; 318 value += scalar * x[iRowP]; 319 y[iColumn] = value; 320 } 321 } else { 322 // skip negative rows 323 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 324 double value = y[iColumn]; 325 int iRowM = indices_[j]; 326 int iRowP = indices_[j + 1]; 327 if (iRowM >= 0) 328 value = scalar * x[iRowM]; 329 if (iRowP >= 0) 330 value += scalar * x[iRowP]; 331 y[iColumn] = value; 332 } 333 } 334 } 335 void ClpNetworkMatrix::times(double scalar, 336 const double *x, double *y, 337 const double * /*rowScale*/, 338 const double * /*columnScale*/) const 339 { 340 // we know it is not scaled 341 times(scalar, x, y); 342 } 343 void ClpNetworkMatrix::transposeTimes(double scalar, 344 const double *x, double *y, 345 const double * /*rowScale*/, 346 const double * /*columnScale*/, 347 double * /*spare*/) const 348 { 349 // we know it is not scaled 350 transposeTimes(scalar, x, y); 355 351 } 356 352 /* Return <code>x * A + y</code> in <code>z</code>. 357 353 Squashes small elements and knows about ClpSimplex */ 358 void 359 ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar, 360 const CoinIndexedVector * rowArray, 361 CoinIndexedVector * y, 362 CoinIndexedVector * columnArray) const 363 { 364 // we know it is not scaled 365 columnArray>clear(); 366 double * pi = rowArray>denseVector(); 367 int numberNonZero = 0; 368 int * index = columnArray>getIndices(); 369 double * array = columnArray>denseVector(); 370 int numberInRowArray = rowArray>getNumElements(); 371 // maybe I need one in OsiSimplex 372 double zeroTolerance = model>zeroTolerance(); 373 int numberRows = model>numberRows(); 354 void ClpNetworkMatrix::transposeTimes(const ClpSimplex *model, double scalar, 355 const CoinIndexedVector *rowArray, 356 CoinIndexedVector *y, 357 CoinIndexedVector *columnArray) const 358 { 359 // we know it is not scaled 360 columnArray>clear(); 361 double *pi = rowArray>denseVector(); 362 int numberNonZero = 0; 363 int *index = columnArray>getIndices(); 364 double *array = columnArray>denseVector(); 365 int numberInRowArray = rowArray>getNumElements(); 366 // maybe I need one in OsiSimplex 367 double zeroTolerance = model>zeroTolerance(); 368 int numberRows = model>numberRows(); 374 369 #ifndef NO_RTTI 375 ClpPlusMinusOneMatrix* rowCopy = 376 dynamic_cast< ClpPlusMinusOneMatrix*>(model>rowCopy()); 370 ClpPlusMinusOneMatrix *rowCopy = dynamic_cast< ClpPlusMinusOneMatrix * >(model>rowCopy()); 377 371 #else 378 ClpPlusMinusOneMatrix* rowCopy = 379 static_cast< ClpPlusMinusOneMatrix*>(model>rowCopy()); 372 ClpPlusMinusOneMatrix *rowCopy = static_cast< ClpPlusMinusOneMatrix * >(model>rowCopy()); 380 373 #endif 381 bool packed = rowArray>packedMode(); 382 double factor = 0.3; 383 // We may not want to do by row if there may be cache problems 384 int numberColumns = model>numberColumns(); 385 // It would be nice to find L2 cache size  for moment 512K 386 // Be slightly optimistic 387 if (numberColumns * sizeof(double) > 1000000) { 388 if (numberRows * 10 < numberColumns) 389 factor = 0.1; 390 else if (numberRows * 4 < numberColumns) 391 factor = 0.15; 392 else if (numberRows * 2 < numberColumns) 393 factor = 0.2; 394 //if (model>numberIterations()%50==0) 395 //printf("%d nonzero\n",numberInRowArray); 396 } 397 if (numberInRowArray > factor * numberRows  !rowCopy) { 398 // do by column 399 int iColumn; 400 assert (!y>getNumElements()); 401 CoinBigIndex j = 0; 402 if (packed) { 403 // need to expand pi into y 404 assert(y>capacity() >= numberRows); 405 double * piOld = pi; 406 pi = y>denseVector(); 407 const int * whichRow = rowArray>getIndices(); 408 int i; 409 // modify pi so can collapse to one loop 410 for (i = 0; i < numberInRowArray; i++) { 411 int iRow = whichRow[i]; 412 pi[iRow] = scalar * piOld[i]; 413 } 414 if (trueNetwork_) { 415 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 416 double value = 0.0; 417 int iRowM = indices_[j]; 418 int iRowP = indices_[j+1]; 419 value = pi[iRowM]; 420 value += pi[iRowP]; 421 if (fabs(value) > zeroTolerance) { 422 array[numberNonZero] = value; 423 index[numberNonZero++] = iColumn; 424 } 425 } 426 } else { 427 // skip negative rows 428 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 429 double value = 0.0; 430 int iRowM = indices_[j]; 431 int iRowP = indices_[j+1]; 432 if (iRowM >= 0) 433 value = pi[iRowM]; 434 if (iRowP >= 0) 435 value += pi[iRowP]; 436 if (fabs(value) > zeroTolerance) { 437 array[numberNonZero] = value; 438 index[numberNonZero++] = iColumn; 439 } 440 } 441 } 442 for (i = 0; i < numberInRowArray; i++) { 443 int iRow = whichRow[i]; 444 pi[iRow] = 0.0; 445 } 446 } else { 447 if (trueNetwork_) { 448 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 449 double value = 0.0; 450 int iRowM = indices_[j]; 451 int iRowP = indices_[j+1]; 452 value = scalar * pi[iRowM]; 453 value += scalar * pi[iRowP]; 454 if (fabs(value) > zeroTolerance) { 455 index[numberNonZero++] = iColumn; 456 array[iColumn] = value; 457 } 458 } 459 } else { 460 // skip negative rows 461 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 462 double value = 0.0; 463 int iRowM = indices_[j]; 464 int iRowP = indices_[j+1]; 465 if (iRowM >= 0) 466 value = scalar * pi[iRowM]; 467 if (iRowP >= 0) 468 value += scalar * pi[iRowP]; 469 if (fabs(value) > zeroTolerance) { 470 index[numberNonZero++] = iColumn; 471 array[iColumn] = value; 472 } 473 } 474 } 374 bool packed = rowArray>packedMode(); 375 double factor = 0.3; 376 // We may not want to do by row if there may be cache problems 377 int numberColumns = model>numberColumns(); 378 // It would be nice to find L2 cache size  for moment 512K 379 // Be slightly optimistic 380 if (numberColumns * sizeof(double) > 1000000) { 381 if (numberRows * 10 < numberColumns) 382 factor = 0.1; 383 else if (numberRows * 4 < numberColumns) 384 factor = 0.15; 385 else if (numberRows * 2 < numberColumns) 386 factor = 0.2; 387 //if (model>numberIterations()%50==0) 388 //printf("%d nonzero\n",numberInRowArray); 389 } 390 if (numberInRowArray > factor * numberRows  !rowCopy) { 391 // do by column 392 int iColumn; 393 assert(!y>getNumElements()); 394 CoinBigIndex j = 0; 395 if (packed) { 396 // need to expand pi into y 397 assert(y>capacity() >= numberRows); 398 double *piOld = pi; 399 pi = y>denseVector(); 400 const int *whichRow = rowArray>getIndices(); 401 int i; 402 // modify pi so can collapse to one loop 403 for (i = 0; i < numberInRowArray; i++) { 404 int iRow = whichRow[i]; 405 pi[iRow] = scalar * piOld[i]; 406 } 407 if (trueNetwork_) { 408 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 409 double value = 0.0; 410 int iRowM = indices_[j]; 411 int iRowP = indices_[j + 1]; 412 value = pi[iRowM]; 413 value += pi[iRowP]; 414 if (fabs(value) > zeroTolerance) { 415 array[numberNonZero] = value; 416 index[numberNonZero++] = iColumn; 475 417 } 476 columnArray>setNumElements(numberNonZero); 477 } else { 478 // do by row 479 rowCopy>transposeTimesByRow(model, scalar, rowArray, y, columnArray); 480 } 418 } 419 } else { 420 // skip negative rows 421 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 422 double value = 0.0; 423 int iRowM = indices_[j]; 424 int iRowP = indices_[j + 1]; 425 if (iRowM >= 0) 426 value = pi[iRowM]; 427 if (iRowP >= 0) 428 value += pi[iRowP]; 429 if (fabs(value) > zeroTolerance) { 430 array[numberNonZero] = value; 431 index[numberNonZero++] = iColumn; 432 } 433 } 434 } 435 for (i = 0; i < numberInRowArray; i++) { 436 int iRow = whichRow[i]; 437 pi[iRow] = 0.0; 438 } 439 } else { 440 if (trueNetwork_) { 441 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 442 double value = 0.0; 443 int iRowM = indices_[j]; 444 int iRowP = indices_[j + 1]; 445 value = scalar * pi[iRowM]; 446 value += scalar * pi[iRowP]; 447 if (fabs(value) > zeroTolerance) { 448 index[numberNonZero++] = iColumn; 449 array[iColumn] = value; 450 } 451 } 452 } else { 453 // skip negative rows 454 for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) { 455 double value = 0.0; 456 int iRowM = indices_[j]; 457 int iRowP = indices_[j + 1]; 458 if (iRowM >= 0) 459 value = scalar * pi[iRowM]; 460 if (iRowP >= 0) 461 value += scalar * pi[iRowP]; 462 if (fabs(value) > zeroTolerance) { 463 index[numberNonZero++] = iColumn; 464 array[iColumn] = value; 465 } 466 } 467 } 468 } 469 columnArray>setNumElements(numberNonZero); 470 } else { 471 // do by row 472 rowCopy>transposeTimesByRow(model, scalar, rowArray, y, columnArray); 473 } 481 474 } 482 475 /* Return <code>x *A in <code>z</code> but 483 476 just for indices in y. */ 484 void 485 ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/, 486 const CoinIndexedVector * rowArray, 487 const CoinIndexedVector * y, 488 CoinIndexedVector * columnArray) const 489 { 490 columnArray>clear(); 491 double * pi = rowArray>denseVector(); 492 double * array = columnArray>denseVector(); 493 int jColumn; 494 int numberToDo = y>getNumElements(); 495 const int * which = y>getIndices(); 496 assert (!rowArray>packedMode()); 497 columnArray>setPacked(); 498 if (trueNetwork_) { 499 for (jColumn = 0; jColumn < numberToDo; jColumn++) { 500 int iColumn = which[jColumn]; 501 double value = 0.0; 502 CoinBigIndex j = iColumn << 1; 503 int iRowM = indices_[j]; 504 int iRowP = indices_[j+1]; 505 value = pi[iRowM]; 506 value += pi[iRowP]; 507 array[jColumn] = value; 508 } 509 } else { 510 // skip negative rows 511 for (jColumn = 0; jColumn < numberToDo; jColumn++) { 512 int iColumn = which[jColumn]; 513 double value = 0.0; 514 CoinBigIndex j = iColumn << 1; 515 int iRowM = indices_[j]; 516 int iRowP = indices_[j+1]; 517 if (iRowM >= 0) 518 value = pi[iRowM]; 519 if (iRowP >= 0) 520 value += pi[iRowP]; 521 array[jColumn] = value; 522 } 523 } 477 void ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/, 478 const CoinIndexedVector *rowArray, 479 const CoinIndexedVector *y, 480 CoinIndexedVector *columnArray) const 481 { 482 columnArray>clear(); 483 double *pi = rowArray>denseVector(); 484 double *array = columnArray>denseVector(); 485 int jColumn; 486 int numberToDo = y>getNumElements(); 487 const int *which = y>getIndices(); 488 assert(!rowArray>packedMode()); 489 columnArray>setPacked(); 490 if (trueNetwork_) { 491 for (jColumn = 0; jColumn < numberToDo; jColumn++) { 492 int iColumn = which[jColumn]; 493 double value = 0.0; 494 CoinBigIndex j = iColumn << 1; 495 int iRowM = indices_[j]; 496 int iRowP = indices_[j + 1]; 497 value = pi[iRowM]; 498 value += pi[iRowP]; 499 array[jColumn] = value; 500 } 501 } else { 502 // skip negative rows 503 for (jColumn = 0; jColumn < numberToDo; jColumn++) { 504 int iColumn = which[jColumn]; 505 double value = 0.0; 506 CoinBigIndex j = iColumn << 1; 507 int iRowM = indices_[j]; 508 int iRowP = indices_[j + 1]; 509 if (iRowM >= 0) 510 value = pi[iRowM]; 511 if (iRowP >= 0) 512 value += pi[iRowP]; 513 array[jColumn] = value; 514 } 515 } 524 516 } 525 517 /// returns number of elements in column part of basis, 526 int 527 ClpNetworkMatrix::countBasis( const int * whichColumn, 528 int & numberColumnBasic) 529 { 530 int i; 531 CoinBigIndex numberElements = 0; 532 if (trueNetwork_) { 533 numberElements = 2 * numberColumnBasic; 534 } else { 535 for (i = 0; i < numberColumnBasic; i++) { 536 int iColumn = whichColumn[i]; 537 CoinBigIndex j = iColumn << 1; 538 int iRowM = indices_[j]; 539 int iRowP = indices_[j+1]; 540 if (iRowM >= 0) 541 numberElements ++; 542 if (iRowP >= 0) 543 numberElements ++; 544 } 545 } 518 int ClpNetworkMatrix::countBasis(const int *whichColumn, 519 int &numberColumnBasic) 520 { 521 int i; 522 CoinBigIndex numberElements = 0; 523 if (trueNetwork_) { 524 numberElements = 2 * numberColumnBasic; 525 } else { 526 for (i = 0; i < numberColumnBasic; i++) { 527 int iColumn = whichColumn[i]; 528 CoinBigIndex j = iColumn << 1; 529 int iRowM = indices_[j]; 530 int iRowP = indices_[j + 1]; 531 if (iRowM >= 0) 532 numberElements++; 533 if (iRowP >= 0) 534 numberElements++; 535 } 536 } 546 537 #if COIN_BIG_INDEX 547 if (numberElements>COIN_INT_MAX) {548 549 550 538 if (numberElements > COIN_INT_MAX) { 539 printf("Factorization too large\n"); 540 abort(); 541 } 551 542 #endif 552 return static_cast<int>(numberElements); 553 } 554 void 555 ClpNetworkMatrix::fillBasis(ClpSimplex * /*model*/, 556 const int * whichColumn, 557 int & numberColumnBasic, 558 int * indexRowU, int * start, 559 int * rowCount, int * columnCount, 560 CoinFactorizationDouble * elementU) 561 { 562 int i; 563 CoinBigIndex numberElements = start[0]; 564 if (trueNetwork_) { 565 for (i = 0; i < numberColumnBasic; i++) { 566 int iColumn = whichColumn[i]; 567 CoinBigIndex j = iColumn << 1; 568 int iRowM = indices_[j]; 569 int iRowP = indices_[j+1]; 570 indexRowU[numberElements] = iRowM; 571 rowCount[iRowM]++; 572 elementU[numberElements] = 1.0; 573 indexRowU[numberElements+1] = iRowP; 574 rowCount[iRowP]++; 575 elementU[numberElements+1] = 1.0; 576 numberElements += 2; 577 start[i+1] = static_cast<int>(numberElements); 578 columnCount[i] = 2; 579 } 580 } else { 581 for (i = 0; i < numberColumnBasic; i++) { 582 int iColumn = whichColumn[i]; 583 CoinBigIndex j = iColumn << 1; 584 int iRowM = indices_[j]; 585 int iRowP = indices_[j+1]; 586 if (iRowM >= 0) { 587 indexRowU[numberElements] = iRowM; 588 rowCount[iRowM]++; 589 elementU[numberElements++] = 1.0; 590 } 591 if (iRowP >= 0) { 592 indexRowU[numberElements] = iRowP; 593 rowCount[iRowP]++; 594 elementU[numberElements++] = 1.0; 595 } 596 start[i+1] = static_cast<int>(numberElements); 597 columnCount[i] = static_cast<int>(numberElements)  start[i]; 598 } 599 } 543 return static_cast< int >(numberElements); 544 } 545 void ClpNetworkMatrix::fillBasis(ClpSimplex * /*model*/, 546 const int *whichColumn, 547 int &numberColumnBasic, 548 int *indexRowU, int *start, 549 int *rowCount, int *columnCount, 550 CoinFactorizationDouble *elementU) 551 { 552 int i; 553 CoinBigIndex numberElements = start[0]; 554 if (trueNetwork_) { 555 for (i = 0; i < numberColumnBasic; i++) { 556 int iColumn = whichColumn[i]; 557 CoinBigIndex j = iColumn << 1; 558 int iRowM = indices_[j]; 559 int iRowP = indices_[j + 1]; 560 indexRowU[numberElements] = iRowM; 561 rowCount[iRowM]++; 562 elementU[numberElements] = 1.0; 563 indexRowU[numberElements + 1] = iRowP; 564 rowCount[iRowP]++; 565 elementU[numberElements + 1] = 1.0; 566 numberElements += 2; 567 start[i + 1] = static_cast< int >(numberElements); 568 columnCount[i] = 2; 569 } 570 } else { 571 for (i = 0; i < numberColumnBasic; i++) { 572 int iColumn = whichColumn[i]; 573 CoinBigIndex j = iColumn << 1; 574 int iRowM = indices_[j]; 575 int iRowP = indices_[j + 1]; 576 if (iRowM >= 0) { 577 indexRowU[numberElements] = iRowM; 578 rowCount[iRowM]++; 579 elementU[numberElements++] = 1.0; 580 } 581 if (iRowP >= 0) { 582 indexRowU[numberElements] = iRowP; 583 rowCount[iRowP]++; 584 elementU[numberElements++] = 1.0; 585 } 586 start[i + 1] = static_cast< int >(numberElements); 587 columnCount[i] = static_cast< int >(numberElements)  start[i]; 588 } 589 } 600 590 #if COIN_BIG_INDEX 601 if (numberElements>COIN_INT_MAX)602 591 if (numberElements > COIN_INT_MAX) 592 abort(); 603 593 #endif 604 594 } 605 595 /* Unpacks a column into an CoinIndexedvector 606 596 */ 607 void 608 ClpNetworkMatrix::unpack(const ClpSimplex * /*model*/, CoinIndexedVector * rowArray, 609 int iColumn) const 610 { 611 CoinBigIndex j = iColumn << 1; 612 int iRowM = indices_[j]; 613 int iRowP = indices_[j+1]; 614 if (iRowM >= 0) 615 rowArray>add(iRowM, 1.0); 616 if (iRowP >= 0) 617 rowArray>add(iRowP, 1.0); 597 void ClpNetworkMatrix::unpack(const ClpSimplex * /*model*/, CoinIndexedVector *rowArray, 598 int iColumn) const 599 { 600 CoinBigIndex j = iColumn << 1; 601 int iRowM = indices_[j]; 602 int iRowP = indices_[j + 1]; 603 if (iRowM >= 0) 604 rowArray>add(iRowM, 1.0); 605 if (iRowP >= 0) 606 rowArray>add(iRowP, 1.0); 618 607 } 619 608 /* Unpacks a column into an CoinIndexedvector … … 621 610 Note that model is NOT const. Bounds and objective could 622 611 be modified if doing column generation (just for this variable) */ 623 void 624 ClpNetworkMatrix::unpackPacked(ClpSimplex * /*model*/, 625 CoinIndexedVector * rowArray, 626 int iColumn) const 627 { 628 int * index = rowArray>getIndices(); 629 double * array = rowArray>denseVector(); 630 int number = 0; 631 CoinBigIndex j = iColumn << 1; 632 int iRowM = indices_[j]; 633 int iRowP = indices_[j+1]; 634 if (iRowM >= 0) { 635 array[number] = 1.0; 636 index[number++] = iRowM; 637 } 638 if (iRowP >= 0) { 639 array[number] = 1.0; 640 index[number++] = iRowP; 641 } 642 rowArray>setNumElements(number); 643 rowArray>setPackedMode(true); 612 void ClpNetworkMatrix::unpackPacked(ClpSimplex * /*model*/, 613 CoinIndexedVector *rowArray, 614 int iColumn) const 615 { 616 int *index = rowArray>getIndices(); 617 double *array = rowArray>denseVector(); 618 int number = 0; 619 CoinBigIndex j = iColumn << 1; 620 int iRowM = indices_[j]; 621 int iRowP = indices_[j + 1]; 622 if (iRowM >= 0) { 623 array[number] = 1.0; 624 index[number++] = iRowM; 625 } 626 if (iRowP >= 0) { 627 array[number] = 1.0; 628 index[number++] = iRowP; 629 } 630 rowArray>setNumElements(number); 631 rowArray>setPackedMode(true); 644 632 } 645 633 /* Adds multiple of a column into an CoinIndexedvector 646 634 You can use quickAdd to add to vector */ 647 void 648 ClpNetworkMatrix::add(const ClpSimplex * /*model*/, CoinIndexedVector * rowArray, 649 int iColumn, double multiplier) const 650 { 651 CoinBigIndex j = iColumn << 1; 652 int iRowM = indices_[j]; 653 int iRowP = indices_[j+1]; 654 if (iRowM >= 0) 655 rowArray>quickAdd(iRowM, multiplier); 656 if (iRowP >= 0) 657 rowArray>quickAdd(iRowP, multiplier); 635 void ClpNetworkMatrix::add(const ClpSimplex * /*model*/, CoinIndexedVector *rowArray, 636 int iColumn, double multiplier) const 637 { 638 CoinBigIndex j = iColumn << 1; 639 int iRowM = indices_[j]; 640 int iRowP = indices_[j + 1]; 641 if (iRowM >= 0) 642 rowArray>quickAdd(iRowM, multiplier); 643 if (iRowP >= 0) 644 rowArray>quickAdd(iRowP, multiplier); 658 645 } 659 646 /* Adds multiple of a column into an array */ 660 void 661 ClpNetworkMatrix::add(const ClpSimplex * /*model*/, double * array, 662 int iColumn, double multiplier) const 663 { 664 CoinBigIndex j = iColumn << 1; 665 int iRowM = indices_[j]; 666 int iRowP = indices_[j+1]; 667 if (iRowM >= 0) 668 array[iRowM] = multiplier; 669 if (iRowP >= 0) 670 array[iRowP] += multiplier; 647 void ClpNetworkMatrix::add(const ClpSimplex * /*model*/, double *array, 648 int iColumn, double multiplier) const 649 { 650 CoinBigIndex j = iColumn << 1; 651 int iRowM = indices_[j]; 652 int iRowP = indices_[j + 1]; 653 if (iRowM >= 0) 654 array[iRowM] = multiplier; 655 if (iRowP >= 0) 656 array[iRowP] += multiplier; 671 657 } 672 658 … … 675 661 ClpNetworkMatrix::getPackedMatrix() const 676 662 { 677 678 assert(trueNetwork_); // fix later679 680 double * elements = new double[numberElements];681 682 683 684 elements[i+1] = 1.0;685 686 CoinBigIndex * starts = new CoinBigIndex [numberColumns_+1];687 688 689 690 691 delete[] lengths_;692 693 694 int *indices = CoinCopyOfArray(indices_, 2 * numberColumns_);695 696 697 698 699 700 701 assert(!indices);702 assert(!lengths_);703 704 663 if (!matrix_) { 664 assert(trueNetwork_); // fix later 665 int numberElements = 2 * numberColumns_; 666 double *elements = new double[numberElements]; 667 CoinBigIndex i; 668 for (i = 0; i < 2 * numberColumns_; i += 2) { 669 elements[i] = 1.0; 670 elements[i + 1] = 1.0; 671 } 672 CoinBigIndex *starts = new CoinBigIndex[numberColumns_ + 1]; 673 for (i = 0; i < numberColumns_ + 1; i++) { 674 starts[i] = 2 * i; 675 } 676 // use assignMatrix to save space 677 delete[] lengths_; 678 lengths_ = NULL; 679 matrix_ = new CoinPackedMatrix(); 680 int *indices = CoinCopyOfArray(indices_, 2 * numberColumns_); 681 matrix_>assignMatrix(true, numberRows_, numberColumns_, 682 getNumElements(), 683 elements, indices, 684 starts, lengths_); 685 assert(!elements); 686 assert(!starts); 687 assert(!indices); 688 assert(!lengths_); 689 } 690 return matrix_; 705 691 } 706 692 /* A vector containing the elements in the packed matrix. Note that there … … 711 697 ClpNetworkMatrix::getElements() const 712 698 { 713 714 715 699 if (!matrix_) 700 getPackedMatrix(); 701 return matrix_>getElements(); 716 702 } 717 703 … … 719 705 ClpNetworkMatrix::getVectorStarts() const 720 706 { 721 722 723 707 if (!matrix_) 708 getPackedMatrix(); 709 return matrix_>getVectorStarts(); 724 710 } 725 711 /* The lengths of the majordimension vectors. */ … … 727 713 ClpNetworkMatrix::getVectorLengths() const 728 714 { 729 assert(trueNetwork_); // fix later730 731 lengths_ = new int[numberColumns_];732 733 734 735 736 737 715 assert(trueNetwork_); // fix later 716 if (!lengths_) { 717 lengths_ = new int[numberColumns_]; 718 int i; 719 for (i = 0; i < numberColumns_; i++) { 720 lengths_[i] = 2; 721 } 722 } 723 return lengths_; 738 724 } 739 725 /* Delete the columns whose indices are listed in <code>indDel</code>. */ 740 void ClpNetworkMatrix::deleteCols(const int numDel, const int * 741 { 742 assert(trueNetwork_);743 744 745 746 char *which = new char[numberColumns_];747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 delete[] lengths_;765 766 767 768 769 int * newIndices = new int[newSize];770 771 772 773 774 775 776 777 778 779 780 delete[] which;781 delete[] indices_;782 783 726 void ClpNetworkMatrix::deleteCols(const int numDel, const int *indDel) 727 { 728 assert(trueNetwork_); 729 int iColumn; 730 int numberBad = 0; 731 // Use array to make sure we can have duplicates 732 char *which = new char[numberColumns_]; 733 memset(which, 0, numberColumns_); 734 int nDuplicate = 0; 735 for (iColumn = 0; iColumn < numDel; iColumn++) { 736 int jColumn = indDel[iColumn]; 737 if (jColumn < 0  jColumn >= numberColumns_) { 738 numberBad++; 739 } else { 740 if (which[jColumn]) 741 nDuplicate++; 742 else 743 which[jColumn] = 1; 744 } 745 } 746 if (numberBad) 747 throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix"); 748 int newNumber = numberColumns_  numDel + nDuplicate; 749 // Get rid of temporary arrays 750 delete[] lengths_; 751 lengths_ = NULL; 752 delete matrix_; 753 matrix_ = NULL; 754 int newSize = 2 * newNumber; 755 int *newIndices = new int[newSize]; 756 newSize = 0; 757 for (iColumn = 0; iColumn < numberColumns_; iColumn++) { 758 if (!which[iColumn]) { 759 CoinBigIndex start; 760 CoinBigIndex i; 761 start = 2 * iColumn; 762 for (i = start; i < start + 2; i++) 763 newIndices[newSize++] = indices_[i]; 764 } 765 } 766 delete[] which; 767 delete[] indices_; 768 indices_ = newIndices; 769 numberColumns_ = newNumber; 784 770 } 785 771 /* Delete the rows whose indices are listed in <code>indDel</code>. */ 786 void ClpNetworkMatrix::deleteRows(const int numDel, const int * 787 { 788 789 790 791 int * which = new int[numberRows_];792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 delete[] which;834 772 void ClpNetworkMatrix::deleteRows(const int numDel, const int *indDel) 773 { 774 int iRow; 775 int numberBad = 0; 776 // Use array to make sure we can have duplicates 777 int *which = new int[numberRows_]; 778 memset(which, 0, numberRows_ * sizeof(int)); 779 for (iRow = 0; iRow < numDel; iRow++) { 780 int jRow = indDel[iRow]; 781 if (jRow < 0  jRow >= numberRows_) { 782 numberBad++; 783 } else { 784 which[jRow] = 1; 785 } 786 } 787 if (numberBad) 788 throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix"); 789 // Only valid of all columns have 0 entries 790 int iColumn; 791 for (iColumn = 0; iColumn < numberColumns_; iColumn++) { 792 CoinBigIndex start; 793 CoinBigIndex i; 794 start = 2 * iColumn; 795 for (i = start; i < start + 2; i++) { 796 int iRow = indices_[i]; 797 if (which[iRow]) 798 numberBad++; 799 } 800 } 801 if (numberBad) 802 throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix"); 803 int newNumber = 0; 804 for (iRow = 0; iRow < numberRows_; iRow++) { 805 if (!which[iRow]) 806 which[iRow] = newNumber++; 807 else 808 which[iRow] = 1; 809 } 810 for (iColumn = 0; iColumn < numberColumns_; iColumn++) { 811 CoinBigIndex start; 812 CoinBigIndex i; 813 start = 2 * iColumn; 814 for (i = start; i < start + 2; i++) { 815 int iRow = indices_[i]; 816 indices_[i] = which[iRow]; 817 } 818 } 819 delete[] which; 820 numberRows_ = newNumber; 835 821 } 836 822 /* Given positive integer weights for each row fills in sum of weights … … 839 825 */ 840 826 CoinBigIndex * 841 ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model, int *inputWeights) const842 { 843 844 845 846 CoinBigIndex *weights = new CoinBigIndex[number];847 848 849 850 851 852 int iRowP = indices_[j+1];853 854 855 856 857 858 859 860 861 862 weights[i+numberColumns] = inputWeights[i];863 864 827 ClpNetworkMatrix::dubiousWeights(const ClpSimplex *model, int *inputWeights) const 828 { 829 int numberRows = model>numberRows(); 830 int numberColumns = model>numberColumns(); 831 int number = numberRows + numberColumns; 832 CoinBigIndex *weights = new CoinBigIndex[number]; 833 int i; 834 for (i = 0; i < numberColumns; i++) { 835 CoinBigIndex j = i << 1; 836 CoinBigIndex count = 0; 837 int iRowM = indices_[j]; 838 int iRowP = indices_[j + 1]; 839 if (iRowM >= 0) { 840 count += inputWeights[iRowM]; 841 } 842 if (iRowP >= 0) { 843 count += inputWeights[iRowP]; 844 } 845 weights[i] = count; 846 } 847 for (i = 0; i < numberRows; i++) { 848 weights[i + numberColumns] = inputWeights[i]; 849 } 850 return weights; 865 851 } 866 852 /* Returns largest and smallest elements of both signs. 867 853 Largest refers to largest absolute value. 868 854 */ 869 void 870 ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative, 871 double & smallestPositive, double & largestPositive) 872 { 873 smallestNegative = 1.0; 874 largestNegative = 1.0; 875 smallestPositive = 1.0; 876 largestPositive = 1.0; 855 void ClpNetworkMatrix::rangeOfElements(double &smallestNegative, double &largestNegative, 856 double &smallestPositive, double &largestPositive) 857 { 858 smallestNegative = 1.0; 859 largestNegative = 1.0; 860 smallestPositive = 1.0; 861 largestPositive = 1.0; 877 862 } 878 863 // Says whether it can do partial pricing 879 bool 880 ClpNetworkMatrix::canDoPartialPricing() const 881 { 882 return true; 864 bool ClpNetworkMatrix::canDoPartialPricing() const 865 { 866 return true; 883 867 } 884 868 // Partial pricing 885 void 886 ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, 887 int & bestSequence, int & numberWanted) 888 { 889 numberWanted = currentWanted_; 890 int j; 891 int start = static_cast<int> (startFraction * numberColumns_); 892 int end = CoinMin(static_cast<int> (endFraction * numberColumns_ + 1), numberColumns_); 893 double tolerance = model>currentDualTolerance(); 894 double * reducedCost = model>djRegion(); 895 const double * duals = model>dualRowSolution(); 896 const double * cost = model>costRegion(); 897 double bestDj; 898 if (bestSequence >= 0) 899 bestDj = fabs(reducedCost[bestSequence]); 900 else 901 bestDj = tolerance; 902 int sequenceOut = model>sequenceOut(); 903 int saveSequence = bestSequence; 904 if (!trueNetwork_) { 905 // Not true network 906 int iSequence; 907 for (iSequence = start; iSequence < end; iSequence++) { 908 if (iSequence != sequenceOut) { 909 double value; 910 int iRowM, iRowP; 911 ClpSimplex::Status status = model>getStatus(iSequence); 912 913 switch(status) { 914 915 case ClpSimplex::basic: 916 case ClpSimplex::isFixed: 917 break; 918 case ClpSimplex::isFree: 919 case ClpSimplex::superBasic: 920 value = cost[iSequence]; 921 j = iSequence << 1; 922 // skip negative rows 923 iRowM = indices_[j]; 924 iRowP = indices_[j+1]; 925 if (iRowM >= 0) 926 value += duals[iRowM]; 927 if (iRowP >= 0) 928 value = duals[iRowP]; 929 value = fabs(value); 930 if (value > FREE_ACCEPT * tolerance) { 931 numberWanted; 932 // we are going to bias towards free (but only if reasonable) 933 value *= FREE_BIAS; 934 if (value > bestDj) { 935 // check flagged variable and correct dj 936 if (!model>flagged(iSequence)) { 937 bestDj = value; 938 bestSequence = iSequence; 939 } else { 940 // just to make sure we don't exit before got something 941 numberWanted++; 942 } 943 } 944 } 945 break; 946 case ClpSimplex::atUpperBound: 947 value = cost[iSequence]; 948 j = iSequence << 1; 949 // skip negative rows 950 iRowM = indices_[j]; 951 iRowP = indices_[j+1]; 952 if (iRowM >= 0) 953 value += duals[iRowM]; 954 if (iRowP >= 0) 955 value = duals[iRowP]; 956 if (value > tolerance) { 957 numberWanted; 958 if (value > bestDj) { 959 // check flagged variable and correct dj 960 if (!model>flagged(iSequence)) { 961 bestDj = value; 962 bestSequence = iSequence; 963 } else { 964 // just to make sure we don't exit before got something 965 numberWanted++; 966 } 967 } 968 } 969 break; 970 case ClpSimplex::atLowerBound: 971 value = cost[iSequence]; 972 j = iSequence << 1; 973 // skip negative rows 974 iRowM = indices_[j]; 975 iRowP = indices_[j+1]; 976 if (iRowM >= 0) 977 value += duals[iRowM]; 978 if (iRowP >= 0) 979 value = duals[iRowP]; 980 value = value; 981 if (value > tolerance) { 982 numberWanted; 983 if (value > bestDj) { 984 // check flagged variable and correct dj 985 if (!model>flagged(iSequence)) { 986 bestDj = value; 987 bestSequence = iSequence; 988 } else { 989 // just to make sure we don't exit before got something 990 numberWanted++; 991 } 992 } 993 } 994 break; 995 } 996 } 997 if (!numberWanted) 998 break; 869 void ClpNetworkMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction, 870 int &bestSequence, int &numberWanted) 871 { 872 numberWanted = currentWanted_; 873 int j; 874 int start = static_cast< int >(startFraction * numberColumns_); 875 int end = CoinMin(static_cast< int >(endFraction * numberColumns_ + 1), numberColumns_); 876 double tolerance = model>currentDualTolerance(); 877 double *reducedCost = model>djRegion(); 878 const double *duals = model>dualRowSolution(); 879 const double *cost = model>costRegion(); 880 double bestDj; 881 if (bestSequence >= 0) 882 bestDj = fabs(reducedCost[bestSequence]); 883 else 884 bestDj = tolerance; 885 int sequenceOut = model>sequenceOut(); 886 int saveSequence = bestSequence; 887 if (!trueNetwork_) { 888 // Not true network 889 int iSequence; 890 for (iSequence = start; iSequence < end; iSequence++) { 891 if (iSequence != sequenceOut) { 892 double value; 893 int iRowM, iRowP; 894 ClpSimplex::Status status = model>getStatus(iSequence); 895 896 switch (status) { 897 898 case ClpSimplex::basic: 899 case ClpSimplex::isFixed: 900 break; 901 case ClpSimplex::isFree: 902 case ClpSimplex::superBasic: 903 value = cost[iSequence]; 904 j = iSequence << 1; 905 // skip negative rows 906 iRowM = indices_[j]; 907 iRowP = indices_[j + 1]; 908 if (iRowM >= 0) 909 value += duals[iRowM]; 910 if (iRowP >= 0) 911 value = duals[iRowP]; 912 value = fabs(value); 913 if (value > FREE_ACCEPT * tolerance) { 914 numberWanted; 915 // we are going to bias towards free (but only if reasonable) 916 value *= FREE_BIAS; 917 if (value > bestDj) { 918 // check flagged variable and correct dj 919 if (!model>flagged(iSequence)) { 920 bestDj = value; 921 bestSequence = iSequence; 922 } else { 923 // just to make sure we don't exit before got something 924 numberWanted++; 925 } 926 } 999 927 } 1000 if (bestSequence != saveSequence) { 1001 // recompute dj 1002 double value = cost[bestSequence]; 1003 j = bestSequence << 1; 1004 // skip negative rows 1005 int iRowM = indices_[j]; 1006 int iRowP = indices_[j+1]; 1007 if (iRowM >= 0) 1008 value += duals[iRowM]; 1009 if (iRowP >= 0) 1010 value = duals[iRowP]; 1011 reducedCost[bestSequence] = value; 1012 savedBestSequence_ = bestSequence; 1013 savedBestDj_ = reducedCost[savedBestSequence_]; 928 break; 929 case ClpSimplex::atUpperBound: 930 value = cost[iSequence]; 931 j = iSequence << 1; 932 // skip negative rows 933 iRowM = indices_[j]; 934 iRowP = indices_[j + 1]; 935 if (iRowM >= 0) 936 value += duals[iRowM]; 937 if (iRowP >= 0) 938 value = duals[iRowP]; 939 if (value > tolerance) { 940 numberWanted; 941 if (value > bestDj) { 942 // check flagged variable and correct dj 943 if (!model>flagged(iSequence)) { 944 bestDj = value; 945 bestSequence = iSequence; 946 } else { 947 // just to make sure we don't exit before got something 948 numberWanted++; 949 } 950 } 1014 951 } 1015 } else { 1016 // true network 1017 int iSequence; 1018 for (iSequence = start; iSequence < end; iSequence++) { 1019 if (iSequence != sequenceOut) { 1020 double value; 1021 int iRowM, iRowP; 1022 ClpSimplex::Status status = model>getStatus(iSequence); 1023 1024 switch(status) { 1025 1026 case ClpSimplex::basic: 1027 case ClpSimplex::isFixed: 1028 break; 1029 case ClpSimplex::isFree: 1030 case ClpSimplex::superBasic: 1031 value = cost[iSequence]; 1032 j = iSequence << 1; 1033 iRowM = indices_[j]; 1034 iRowP = indices_[j+1]; 1035 value += duals[iRowM]; 1036 value = duals[iRowP]; 1037 value = fabs(value); 1038 if (value > FREE_ACCEPT * tolerance) { 1039 numberWanted; 1040 // we are going to bias towards free (but only if reasonable) 1041 value *= FREE_BIAS; 1042 if (value > bestDj) { 1043 // check flagged variable and correct dj 1044 if (!model>flagged(iSequence)) { 1045 bestDj = value; 1046 bestSequence = iSequence; 1047 } else { 1048 // just to make sure we don't exit before got something 1049 numberWanted++; 1050 } 1051 } 1052 } 1053 break; 1054 case ClpSimplex::atUpperBound: 1055 value = cost[iSequence]; 1056 j = iSequence << 1; 1057 iRowM = indices_[j]; 1058 iRowP = indices_[j+1]; 1059 value += duals[iRowM]; 1060 value = duals[iRowP]; 1061 if (value > tolerance) { 1062 numberWanted; 1063 if (value > bestDj) { 1064 // check flagged variable and correct dj 1065 if (!model>flagged(iSequence)) { 1066 bestDj = value; 1067 bestSequence = iSequence; 1068 } else { 1069 // just to make sure we don't exit before got something 1070 numberWanted++; 1071 } 1072 } 1073 } 1074 break; 1075 case ClpSimplex::atLowerBound: 1076 value = cost[iSequence]; 1077 j = iSequence << 1; 1078 iRowM = indices_[j]; 1079 iRowP = indices_[j+1]; 1080 value += duals[iRowM]; 1081 value = duals[iRowP]; 1082 value = value; 1083 if (value > tolerance) { 1084 numberWanted; 1085 if (value > bestDj) { 1086 // check flagged variable and correct dj 1087 if (!model>flagged(iSequence)) { 1088 bestDj = value; 1089 bestSequence = iSequence; 1090 } else { 1091 // just to make sure we don't exit before got something 1092 numberWanted++; 1093 } 1094 } 1095 } 1096 break; 1097 } 1098 } 1099 if (!numberWanted) 1100 break; 952 break; 953 case ClpSimplex::atLowerBound: 954 value = cost[iSequence]; 955 j = iSequence << 1; 956 // skip negative rows 957 iRowM = indices_[j]; 958 iRowP = indices_[j + 1]; 959 if (iRowM >= 0) 960 value += duals[iRowM]; 961 if (iRowP >= 0) 962 value = duals[iRowP]; 963 value = value; 964 if (value > tolerance) { 965 numberWanted; 966 if (value > bestDj) { 967 // check flagged variable and correct dj 968 if (!model>flagged(iSequence)) { 969 bestDj = value; 970 bestSequence = iSequence; 971 } else { 972 // just to make sure we don't exit before got something 973 numberWanted++; 974 } 975 } 1101 976 } 1102 if (bestSequence != saveSequence) { 1103 // recompute dj 1104 double value = cost[bestSequence]; 1105 j = bestSequence << 1; 1106 int iRowM = indices_[j]; 1107 int iRowP = indices_[j+1]; 1108 value += duals[iRowM]; 1109 value = duals[iRowP]; 1110 reducedCost[bestSequence] = value; 1111 savedBestSequence_ = bestSequence; 1112 savedBestDj_ = reducedCost[savedBestSequence_]; 977 break; 978 } 979 } 980 if (!numberWanted) 981 break; 982 } 983 if (bestSequence != saveSequence) { 984 // recompute dj 985 double value = cost[bestSequence]; 986 j = bestSequence << 1; 987 // skip negative rows 988 int iRowM = indices_[j]; 989 int iRowP = indices_[j + 1]; 990 if (iRowM >= 0) 991 value += duals[iRowM]; 992 if (iRowP >= 0) 993 value = duals[iRowP]; 994 reducedCost[bestSequence] = value; 995 savedBestSequence_ = bestSequence; 996 savedBestDj_ = reducedCost[savedBestSequence_]; 997 } 998 } else { 999 // true network 1000 int iSequence; 1001 for (iSequence = start; iSequence < end; iSequence++) { 1002 if (iSequence != sequenceOut) { 1003 double value; 1004 int iRowM, iRowP; 1005 ClpSimplex::Status status = model>getStatus(iSequence); 1006 1007 switch (status) { 1008 1009 case ClpSimplex::basic: 1010 case ClpSimplex::isFixed: 1011 break; 1012 case ClpSimplex::isFree: 1013 case ClpSimplex::superBasic: 1014 value = cost[iSequence]; 1015 j = iSequence << 1; 1016 iRowM = indices_[j]; 1017 iRowP = indices_[j + 1]; 1018 value += duals[iRowM]; 1019 value = duals[iRowP]; 1020 value = fabs(value); 1021 if (value > FREE_ACCEPT * tolerance) { 1022 numberWanted; 1023 // we are going to bias towards free (but only if reasonable) 1024 value *= FREE_BIAS; 1025 if (value > bestDj) { 1026 // check flagged variable and correct dj 1027 if (!model>flagged(iSequence)) { 1028 bestDj = value; 1029 bestSequence = iSequence; 1030 } else { 1031 // just to make sure we don't exit before got something 1032 numberWanted++; 1033 } 1034 } 1113 1035 } 1114 } 1115 currentWanted_ = numberWanted; 1036 break; 1037 case ClpSimplex::atUpperBound: 1038 value = cost[iSequence]; 1039 j = iSequence << 1; 1040 iRowM = indices_[j]; 1041 iRowP = indices_[j + 1]; 1042 value += duals[iRowM]; 1043 value = duals[iRowP]; 1044 if (value > tolerance) { 1045 numberWanted; 1046 if (value > bestDj) { 1047 // check flagged variable and correct dj 1048 if (!model>flagged(iSequence)) { 1049 bestDj = value; 1050 bestSequence = iSequence; 1051 } else { 1052 // just to make sure we don't exit before got something 1053 numberWanted++; 1054 } 1055 } 1056 } 1057 break; 1058 case ClpSimplex::atLowerBound: 1059 value = cost[iSequence]; 1060 j = iSequence << 1; 1061 iRowM = indices_[j]; 1062 iRowP = indices_[j + 1]; 1063 value += duals[iRowM]; 1064 value = duals[iRowP]; 1065 value = value; 1066 if (value > tolerance) { 1067 numberWanted; 1068 if (value > bestDj) { 1069 // check flagged variable and correct dj 1070 if (!model>flagged(iSequence)) { 1071 bestDj = value; 1072 bestSequence = iSequence; 1073 } else { 1074 // just to make sure we don't exit before got something 1075 numberWanted++; 1076 } 1077 } 1078 } 1079 break; 1080 } 1081 } 1082 if (!numberWanted) 1083 break; 1084 } 1085 if (bestSequence != saveSequence) { 1086 // recompute dj 1087 double value = cost[bestSequence]; 1088 j = bestSequence << 1; 1089 int iRowM = indices_[j]; 1090 int iRowP = indices_[j + 1]; 1091 value += duals[iRowM]; 1092 value = duals[iRowP]; 1093 reducedCost[bestSequence] = value; 1094 savedBestSequence_ = bestSequence; 1095 savedBestDj_ = reducedCost[savedBestSequence_]; 1096 } 1097 } 1098 currentWanted_ = numberWanted; 1116 1099 } 1117 1100 // Allow any parts of a created CoinMatrix to be deleted 1118 void 1119 ClpNetworkMatrix::releasePackedMatrix() const 1120 { 1121 delete matrix_; 1122 delete [] lengths_; 1123 matrix_ = NULL; 1124 lengths_ = NULL; 1101 void ClpNetworkMatrix::releasePackedMatrix() const 1102 { 1103 delete matrix_; 1104 delete[] lengths_; 1105 matrix_ = NULL; 1106 lengths_ = NULL; 1125 1107 } 1126 1108 // Append Columns 1127 void 1128 ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns) 1129 { 1130 int iColumn; 1131 int numberBad = 0; 1132 for (iColumn = 0; iColumn < number; iColumn++) { 1133 int n = columns[iColumn]>getNumElements(); 1134 const double * element = columns[iColumn]>getElements(); 1135 if (n != 2) 1136 numberBad++; 1137 if (fabs(element[0]) != 1.0  fabs(element[1]) != 1.0) 1138 numberBad++; 1139 else if (element[0]*element[1] != 1.0) 1140 numberBad++; 1141 } 1142 if (numberBad) 1143 throw CoinError("Not network", "appendCols", "ClpNetworkMatrix"); 1144 // Get rid of temporary arrays 1145 delete [] lengths_; 1146 lengths_ = NULL; 1147 delete matrix_; 1148 matrix_ = NULL; 1149 CoinBigIndex size = 2 * number; 1150 int * temp2 = new int [numberColumns_*2+size]; 1151 CoinMemcpyN(indices_, numberColumns_ * 2, temp2); 1152 delete [] indices_; 1153 indices_ = temp2; 1154 // now add 1155 size = 2 * numberColumns_; 1156 for (iColumn = 0; iColumn < number; iColumn++) { 1157 const int * row = columns[iColumn]>getIndices(); 1158 const double * element = columns[iColumn]>getElements(); 1159 if (element[0] == 1.0) { 1160 indices_[size++] = row[0]; 1161 indices_[size++] = row[1]; 1162 } else { 1163 indices_[size++] = row[1]; 1164 indices_[size++] = row[0]; 1165 } 1166 } 1167 1168 numberColumns_ += number; 1109 void ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase *const *columns) 1110 { 1111 int iColumn; 1112 int numberBad = 0; 1113 for (iColumn = 0; iColumn < number; iColumn++) { 1114 int n = columns[iColumn]>getNumElements(); 1115 const double *element = columns[iColumn]>getElements(); 1116 if (n != 2) 1117 numberBad++; 1118 if (fabs(element[0]) != 1.0  fabs(element[1]) != 1.0) 1119 numberBad++; 1120 else if (element[0] * element[1] != 1.0) 1121 numberBad++; 1122 } 1123 if (numberBad) 1124 throw CoinError("Not network", "appendCols", "ClpNetworkMatrix"); 1125 // Get rid of temporary arrays 1126 delete[] lengths_; 1127 lengths_ = NULL; 1128 delete matrix_; 1129 matrix_ = NULL; 1130 CoinBigIndex size = 2 * number; 1131 int *temp2 = new int[numberColumns_ * 2 + size]; 1132 CoinMemcpyN(indices_, numberColumns_ * 2, temp2); 1133 delete[] indices_; 1134 indices_ = temp2; 1135 // now add 1136 size = 2 * numberColumns_; 1137 for (iColumn = 0; iColumn < number; iColumn++) { 1138 const int *row = columns[iColumn]>getIndices(); 1139 const double *element = columns[iColumn]>getElements(); 1140 if (element[0] == 1.0) { 1141 indices_[size++] = row[0]; 1142 indices_[size++] = row[1]; 1143 } else { 1144 indices_[size++] = row[1]; 1145 indices_[size++] = row[0]; 1146 } 1147 } 1148 1149 numberColumns_ += number; 1169 1150 } 1170 1151 // Append Rows 1171 void 1172 ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows) 1173 { 1174 // must be zero arrays 1175 int numberBad = 0; 1176 int iRow; 1177 for (iRow = 0; iRow < number; iRow++) { 1178 numberBad += rows[iRow]>getNumElements(); 1179 } 1180 if (numberBad) 1181 throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix"); 1182 numberRows_ += number; 1152 void ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase *const *rows) 1153 { 1154 // must be zero arrays 1155 int numberBad = 0; 1156 int iRow; 1157 for (iRow = 0; iRow < number; iRow++) { 1158 numberBad += rows[iRow]>getNumElements(); 1159 } 1160 if (numberBad) 1161 throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix"); 1162 numberRows_ += number; 1183 1163 } 1184 1164 #ifndef SLIM_CLP … … 1187 1167 number of columns1/rows1 (if numberOther>0) or duplicates 1188 1168 If 0 then rows, 1 if columns */ 1189 int 1190 ClpNetworkMatrix::appendMatrix(int number, int type, 1191 const CoinBigIndex * starts, const int * index, 1192 const double * element, int /*numberOther*/) 1193 { 1194 int numberErrors = 0; 1195 // make into CoinPackedVector 1196 CoinPackedVectorBase ** vectors = 1197 new CoinPackedVectorBase * [number]; 1198 int iVector; 1199 for (iVector = 0; iVector < number; iVector++) { 1200 CoinBigIndex iStart = starts[iVector]; 1201 vectors[iVector] = 1202 new CoinPackedVector(static_cast<int>(starts[iVector+1]  iStart), 1203 index + iStart, element + iStart); 1204 } 1205 if (type == 0) { 1206 // rows 1207 appendRows(number, vectors); 1208 } else { 1209 // columns 1210 appendCols(number, vectors); 1211 } 1212 for (iVector = 0; iVector < number; iVector++) 1213 delete vectors[iVector]; 1214 delete [] vectors; 1215 return numberErrors; 1169 int ClpNetworkMatrix::appendMatrix(int number, int type, 1170 const CoinBigIndex *starts, const int *index, 1171 const double *element, int /*numberOther*/) 1172 { 1173 int numberErrors = 0; 1174 // make into CoinPackedVector 1175 CoinPackedVectorBase **vectors = new CoinPackedVectorBase *[number]; 1176 int iVector; 1177 for (iVector = 0; iVector < number; iVector++) { 1178 CoinBigIndex iStart = starts[iVector]; 1179 vectors[iVector] = new CoinPackedVector(static_cast< int >(starts[iVector + 1]  iStart), 1180 index + iStart, element + iStart); 1181 } 1182 if (type == 0) { 1183 // rows 1184 appendRows(number, vectors); 1185 } else { 1186 // columns 1187 appendCols(number, vectors); 1188 } 1189 for (iVector = 0; iVector < number; iVector++) 1190 delete vectors[iVector]; 1191 delete[] vectors; 1192 return numberErrors; 1216 1193 } 1217 1194 #endif … … 1219 1196 and order is as given */ 1220 1197 ClpMatrixBase * 1221 ClpNetworkMatrix::subsetClone (int numberRows, const int *whichRows,1222 1223 const int *whichColumns) const1224 { 1225 1226 1198 ClpNetworkMatrix::subsetClone(int numberRows, const int *whichRows, 1199 int numberColumns, 1200 const int *whichColumns) const 1201 { 1202 return new ClpNetworkMatrix(*this, numberRows, whichRows, 1203 numberColumns, whichColumns); 1227 1204 } 1228 1205 /* Subset constructor (without gaps). Duplicates are allowed 1229 1206 and order is as given */ 1230 ClpNetworkMatrix::ClpNetworkMatrix ( 1231 const ClpNetworkMatrix & rhs, 1232 int numberRows, const int * whichRow, 1233 int numberColumns, const int * whichColumn) 1234 : ClpMatrixBase(rhs) 1235 { 1236 setType(11); 1237 matrix_ = NULL; 1238 lengths_ = NULL; 1239 indices_ = new int[2*numberColumns];; 1240 numberRows_ = numberRows; 1241 numberColumns_ = numberColumns; 1242 trueNetwork_ = true; 1243 int iColumn; 1244 int numberBad = 0; 1245 int * which = new int [rhs.numberRows_]; 1246 int iRow; 1247 for (iRow = 0; iRow < rhs.numberRows_; iRow++) 1248 which[iRow] = 1; 1249 int n = 0; 1250 for (iRow = 0; iRow < numberRows; iRow++) { 1251 int jRow = whichRow[iRow]; 1252 assert (jRow >= 0 && jRow < rhs.numberRows_); 1253 which[jRow] = n++; 1254 } 1255 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1256 CoinBigIndex start; 1257 CoinBigIndex i; 1258 start = 2 * iColumn; 1259 CoinBigIndex offset = 2 * whichColumn[iColumn]  start; 1260 for (i = start; i < start + 2; i++) { 1261 int iRow = rhs.indices_[i+offset]; 1262 iRow = which[iRow]; 1263 if (iRow < 0) 1264 numberBad++; 1265 else 1266 indices_[i] = iRow; 1267 } 1268 } 1269 if (numberBad) 1270 throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix"); 1271 } 1207 ClpNetworkMatrix::ClpNetworkMatrix( 1208 const ClpNetworkMatrix &rhs, 1209 int numberRows, const int *whichRow, 1210 int numberColumns, const int *whichColumn) 1211 : ClpMatrixBase(rhs) 1212 { 1213 setType(11); 1214 matrix_ = NULL; 1215 lengths_ = NULL; 1216 indices_ = new int[2 * numberColumns]; 1217 ; 1218 numberRows_ = numberRows; 1219 numberColumns_ = numberColumns; 1220 trueNetwork_ = true; 1221 int iColumn; 1222 int numberBad = 0; 1223 int *which = new int[rhs.numberRows_]; 1224 int iRow; 1225 for (iRow = 0; iRow < rhs.numberRows_; iRow++) 1226 which[iRow] = 1; 1227 int n = 0; 1228 for (iRow = 0; iRow < numberRows; iRow++) { 1229 int jRow = whichRow[iRow]; 1230 assert(jRow >= 0 && jRow < rhs.numberRows_); 1231 which[jRow] = n++; 1232 } 1233 for (iColumn = 0; iColumn < numberColumns; iColumn++) { 1234 CoinBigIndex start; 1235 CoinBigIndex i; 1236 start = 2 * iColumn; 1237 CoinBigIndex offset = 2 * whichColumn[iColumn]  start; 1238 for (i = start; i < start + 2; i++) { 1239 int iRow = rhs.indices_[i + offset]; 1240 iRow = which[iRow]; 1241 if (iRow < 0) 1242 numberBad++; 1243 else 1244 indices_[i] = iRow; 1245 } 1246 } 1247 if (numberBad) 1248 throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix"); 1249 } 1250 1251 /* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2 1252 */
Note: See TracChangeset
for help on using the changeset viewer.