72 | | numberSets_ = rhs.numberSets_; |
73 | | saveNumber_ = rhs.saveNumber_; |
74 | | possiblePivotKey_ = rhs.possiblePivotKey_; |
75 | | gubSlackIn_ = rhs.gubSlackIn_; |
76 | | start_ = ClpCopyOfArray(rhs.start_, numberSets_); |
77 | | end_ = ClpCopyOfArray(rhs.end_, numberSets_); |
78 | | lower_ = ClpCopyOfArray(rhs.lower_, numberSets_); |
79 | | upper_ = ClpCopyOfArray(rhs.upper_, numberSets_); |
80 | | status_ = ClpCopyOfArray(rhs.status_, numberSets_); |
81 | | saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_); |
82 | | savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_); |
83 | | int numberColumns = getNumCols(); |
84 | | backward_ = ClpCopyOfArray(rhs.backward_, numberColumns); |
85 | | backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns); |
86 | | changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_); |
87 | | fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1); |
88 | | keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); |
89 | | // find longest set |
90 | | int * longest = new int[numberSets_]; |
91 | | CoinZeroN(longest, numberSets_); |
92 | | int j; |
93 | | for (j = 0; j < numberColumns; j++) { |
94 | | int iSet = backward_[j]; |
95 | | if (iSet >= 0) |
96 | | longest[iSet]++; |
97 | | } |
98 | | int length = 0; |
99 | | for (j = 0; j < numberSets_; j++) |
100 | | length = CoinMax(length, longest[j]); |
101 | | next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length); |
102 | | toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); |
103 | | sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; |
104 | | sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; |
105 | | sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; |
106 | | sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; |
107 | | infeasibilityWeight_ = rhs.infeasibilityWeight_; |
108 | | numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; |
109 | | numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; |
110 | | noCheck_ = rhs.noCheck_; |
111 | | firstGub_ = rhs.firstGub_; |
112 | | lastGub_ = rhs.lastGub_; |
113 | | gubType_ = rhs.gubType_; |
114 | | model_ = rhs.model_; |
| 71 | numberSets_ = rhs.numberSets_; |
| 72 | saveNumber_ = rhs.saveNumber_; |
| 73 | possiblePivotKey_ = rhs.possiblePivotKey_; |
| 74 | gubSlackIn_ = rhs.gubSlackIn_; |
| 75 | start_ = ClpCopyOfArray(rhs.start_, numberSets_); |
| 76 | end_ = ClpCopyOfArray(rhs.end_, numberSets_); |
| 77 | lower_ = ClpCopyOfArray(rhs.lower_, numberSets_); |
| 78 | upper_ = ClpCopyOfArray(rhs.upper_, numberSets_); |
| 79 | status_ = ClpCopyOfArray(rhs.status_, numberSets_); |
| 80 | saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_); |
| 81 | savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_); |
| 82 | int numberColumns = getNumCols(); |
| 83 | backward_ = ClpCopyOfArray(rhs.backward_, numberColumns); |
| 84 | backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns); |
| 85 | changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_); |
| 86 | fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1); |
| 87 | keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); |
| 88 | // find longest set |
| 89 | int *longest = new int[numberSets_]; |
| 90 | CoinZeroN(longest, numberSets_); |
| 91 | int j; |
| 92 | for (j = 0; j < numberColumns; j++) { |
| 93 | int iSet = backward_[j]; |
| 94 | if (iSet >= 0) |
| 95 | longest[iSet]++; |
| 96 | } |
| 97 | int length = 0; |
| 98 | for (j = 0; j < numberSets_; j++) |
| 99 | length = CoinMax(length, longest[j]); |
| 100 | next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length); |
| 101 | toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); |
| 102 | sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_; |
| 103 | sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; |
| 104 | sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; |
| 105 | sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; |
| 106 | infeasibilityWeight_ = rhs.infeasibilityWeight_; |
| 107 | numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; |
| 108 | numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; |
| 109 | noCheck_ = rhs.noCheck_; |
| 110 | firstGub_ = rhs.firstGub_; |
| 111 | lastGub_ = rhs.lastGub_; |
| 112 | gubType_ = rhs.gubType_; |
| 113 | model_ = rhs.model_; |
192 | | int iSet; |
193 | | for (iSet = 0; iSet < numberSets_; iSet++) { |
194 | | // set key variable as slack |
195 | | keyVariable_[iSet] = iSet + numberColumns; |
196 | | if (start_[iSet] < 0 || start_[iSet] >= numberColumns) |
197 | | throw CoinError("Index out of range", "constructor", "ClpGubMatrix"); |
198 | | if (end_[iSet] < 0 || end_[iSet] > numberColumns) |
199 | | throw CoinError("Index out of range", "constructor", "ClpGubMatrix"); |
200 | | if (end_[iSet] <= start_[iSet]) |
201 | | throw CoinError("Empty or negative set", "constructor", "ClpGubMatrix"); |
202 | | if (start_[iSet] < last) |
203 | | throw CoinError("overlapping or non-monotonic sets", "constructor", "ClpGubMatrix"); |
204 | | last = end_[iSet]; |
205 | | int j; |
206 | | for (j = start_[iSet]; j < end_[iSet]; j++) |
207 | | backward_[j] = iSet; |
208 | | } |
209 | | // Find type of gub |
210 | | firstGub_ = numberColumns + 1; |
211 | | lastGub_ = -1; |
212 | | int i; |
213 | | for (i = 0; i < numberColumns; i++) { |
214 | | if (backward_[i] >= 0) { |
215 | | firstGub_ = CoinMin(firstGub_, i); |
216 | | lastGub_ = CoinMax(lastGub_, i); |
217 | | } |
218 | | } |
219 | | gubType_ = 0; |
220 | | // adjust lastGub_ |
221 | | if (lastGub_ > 0) |
222 | | lastGub_++; |
223 | | for (i = firstGub_; i < lastGub_; i++) { |
224 | | if (backward_[i] < 0) { |
225 | | gubType_ = 1; |
226 | | printf("interior non gub %d\n", i); |
227 | | break; |
228 | | } |
229 | | } |
230 | | if (status) { |
231 | | status_ = ClpCopyOfArray(status, numberSets_); |
232 | | } else { |
233 | | status_ = new unsigned char [numberSets_]; |
234 | | memset(status_, 0, numberSets_); |
235 | | int i; |
236 | | for (i = 0; i < numberSets_; i++) { |
237 | | // make slack key |
238 | | setStatus(i, ClpSimplex::basic); |
239 | | } |
240 | | } |
241 | | saveStatus_ = new unsigned char [numberSets_]; |
242 | | memset(saveStatus_, 0, numberSets_); |
243 | | savedKeyVariable_ = new int [numberSets_]; |
244 | | memset(savedKeyVariable_, 0, numberSets_ * sizeof(int)); |
245 | | noCheck_ = -1; |
246 | | infeasibilityWeight_ = 0.0; |
| 191 | int iSet; |
| 192 | for (iSet = 0; iSet < numberSets_; iSet++) { |
| 193 | // set key variable as slack |
| 194 | keyVariable_[iSet] = iSet + numberColumns; |
| 195 | if (start_[iSet] < 0 || start_[iSet] >= numberColumns) |
| 196 | throw CoinError("Index out of range", "constructor", "ClpGubMatrix"); |
| 197 | if (end_[iSet] < 0 || end_[iSet] > numberColumns) |
| 198 | throw CoinError("Index out of range", "constructor", "ClpGubMatrix"); |
| 199 | if (end_[iSet] <= start_[iSet]) |
| 200 | throw CoinError("Empty or negative set", "constructor", "ClpGubMatrix"); |
| 201 | if (start_[iSet] < last) |
| 202 | throw CoinError("overlapping or non-monotonic sets", "constructor", "ClpGubMatrix"); |
| 203 | last = end_[iSet]; |
| 204 | int j; |
| 205 | for (j = start_[iSet]; j < end_[iSet]; j++) |
| 206 | backward_[j] = iSet; |
| 207 | } |
| 208 | // Find type of gub |
| 209 | firstGub_ = numberColumns + 1; |
| 210 | lastGub_ = -1; |
| 211 | int i; |
| 212 | for (i = 0; i < numberColumns; i++) { |
| 213 | if (backward_[i] >= 0) { |
| 214 | firstGub_ = CoinMin(firstGub_, i); |
| 215 | lastGub_ = CoinMax(lastGub_, i); |
| 216 | } |
| 217 | } |
| 218 | gubType_ = 0; |
| 219 | // adjust lastGub_ |
| 220 | if (lastGub_ > 0) |
| 221 | lastGub_++; |
| 222 | for (i = firstGub_; i < lastGub_; i++) { |
| 223 | if (backward_[i] < 0) { |
| 224 | gubType_ = 1; |
| 225 | printf("interior non gub %d\n", i); |
| 226 | break; |
| 227 | } |
| 228 | } |
| 229 | if (status) { |
| 230 | status_ = ClpCopyOfArray(status, numberSets_); |
| 231 | } else { |
| 232 | status_ = new unsigned char[numberSets_]; |
| 233 | memset(status_, 0, numberSets_); |
| 234 | int i; |
| 235 | for (i = 0; i < numberSets_; i++) { |
| 236 | // make slack key |
| 237 | setStatus(i, ClpSimplex::basic); |
| 238 | } |
| 239 | } |
| 240 | saveStatus_ = new unsigned char[numberSets_]; |
| 241 | memset(saveStatus_, 0, numberSets_); |
| 242 | savedKeyVariable_ = new int[numberSets_]; |
| 243 | memset(savedKeyVariable_, 0, numberSets_ * sizeof(int)); |
| 244 | noCheck_ = -1; |
| 245 | infeasibilityWeight_ = 0.0; |
313 | | if (this != &rhs) { |
314 | | ClpPackedMatrix::operator=(rhs); |
315 | | delete [] start_; |
316 | | delete [] end_; |
317 | | delete [] lower_; |
318 | | delete [] upper_; |
319 | | delete [] status_; |
320 | | delete [] saveStatus_; |
321 | | delete [] savedKeyVariable_; |
322 | | delete [] backward_; |
323 | | delete [] backToPivotRow_; |
324 | | delete [] changeCost_; |
325 | | delete [] keyVariable_; |
326 | | delete [] next_; |
327 | | delete [] toIndex_; |
328 | | delete [] fromIndex_; |
329 | | numberSets_ = rhs.numberSets_; |
330 | | saveNumber_ = rhs.saveNumber_; |
331 | | possiblePivotKey_ = rhs.possiblePivotKey_; |
332 | | gubSlackIn_ = rhs.gubSlackIn_; |
333 | | start_ = ClpCopyOfArray(rhs.start_, numberSets_); |
334 | | end_ = ClpCopyOfArray(rhs.end_, numberSets_); |
335 | | lower_ = ClpCopyOfArray(rhs.lower_, numberSets_); |
336 | | upper_ = ClpCopyOfArray(rhs.upper_, numberSets_); |
337 | | status_ = ClpCopyOfArray(rhs.status_, numberSets_); |
338 | | saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_); |
339 | | savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_); |
340 | | int numberColumns = getNumCols(); |
341 | | backward_ = ClpCopyOfArray(rhs.backward_, numberColumns); |
342 | | backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns); |
343 | | changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_); |
344 | | fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1); |
345 | | keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); |
346 | | // find longest set |
347 | | int * longest = new int[numberSets_]; |
348 | | CoinZeroN(longest, numberSets_); |
349 | | int j; |
350 | | for (j = 0; j < numberColumns; j++) { |
351 | | int iSet = backward_[j]; |
352 | | if (iSet >= 0) |
353 | | longest[iSet]++; |
354 | | } |
355 | | int length = 0; |
356 | | for (j = 0; j < numberSets_; j++) |
357 | | length = CoinMax(length, longest[j]); |
358 | | next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length); |
359 | | toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); |
360 | | sumDualInfeasibilities_ = rhs. sumDualInfeasibilities_; |
361 | | sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; |
362 | | sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; |
363 | | sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; |
364 | | infeasibilityWeight_ = rhs.infeasibilityWeight_; |
365 | | numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; |
366 | | numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; |
367 | | noCheck_ = rhs.noCheck_; |
368 | | firstGub_ = rhs.firstGub_; |
369 | | lastGub_ = rhs.lastGub_; |
370 | | gubType_ = rhs.gubType_; |
371 | | model_ = rhs.model_; |
372 | | } |
373 | | return *this; |
| 311 | if (this != &rhs) { |
| 312 | ClpPackedMatrix::operator=(rhs); |
| 313 | delete[] start_; |
| 314 | delete[] end_; |
| 315 | delete[] lower_; |
| 316 | delete[] upper_; |
| 317 | delete[] status_; |
| 318 | delete[] saveStatus_; |
| 319 | delete[] savedKeyVariable_; |
| 320 | delete[] backward_; |
| 321 | delete[] backToPivotRow_; |
| 322 | delete[] changeCost_; |
| 323 | delete[] keyVariable_; |
| 324 | delete[] next_; |
| 325 | delete[] toIndex_; |
| 326 | delete[] fromIndex_; |
| 327 | numberSets_ = rhs.numberSets_; |
| 328 | saveNumber_ = rhs.saveNumber_; |
| 329 | possiblePivotKey_ = rhs.possiblePivotKey_; |
| 330 | gubSlackIn_ = rhs.gubSlackIn_; |
| 331 | start_ = ClpCopyOfArray(rhs.start_, numberSets_); |
| 332 | end_ = ClpCopyOfArray(rhs.end_, numberSets_); |
| 333 | lower_ = ClpCopyOfArray(rhs.lower_, numberSets_); |
| 334 | upper_ = ClpCopyOfArray(rhs.upper_, numberSets_); |
| 335 | status_ = ClpCopyOfArray(rhs.status_, numberSets_); |
| 336 | saveStatus_ = ClpCopyOfArray(rhs.saveStatus_, numberSets_); |
| 337 | savedKeyVariable_ = ClpCopyOfArray(rhs.savedKeyVariable_, numberSets_); |
| 338 | int numberColumns = getNumCols(); |
| 339 | backward_ = ClpCopyOfArray(rhs.backward_, numberColumns); |
| 340 | backToPivotRow_ = ClpCopyOfArray(rhs.backToPivotRow_, numberColumns); |
| 341 | changeCost_ = ClpCopyOfArray(rhs.changeCost_, getNumRows() + numberSets_); |
| 342 | fromIndex_ = ClpCopyOfArray(rhs.fromIndex_, getNumRows() + numberSets_ + 1); |
| 343 | keyVariable_ = ClpCopyOfArray(rhs.keyVariable_, numberSets_); |
| 344 | // find longest set |
| 345 | int *longest = new int[numberSets_]; |
| 346 | CoinZeroN(longest, numberSets_); |
| 347 | int j; |
| 348 | for (j = 0; j < numberColumns; j++) { |
| 349 | int iSet = backward_[j]; |
| 350 | if (iSet >= 0) |
| 351 | longest[iSet]++; |
| 352 | } |
| 353 | int length = 0; |
| 354 | for (j = 0; j < numberSets_; j++) |
| 355 | length = CoinMax(length, longest[j]); |
| 356 | next_ = ClpCopyOfArray(rhs.next_, numberColumns + numberSets_ + 2 * length); |
| 357 | toIndex_ = ClpCopyOfArray(rhs.toIndex_, numberSets_); |
| 358 | sumDualInfeasibilities_ = rhs.sumDualInfeasibilities_; |
| 359 | sumPrimalInfeasibilities_ = rhs.sumPrimalInfeasibilities_; |
| 360 | sumOfRelaxedDualInfeasibilities_ = rhs.sumOfRelaxedDualInfeasibilities_; |
| 361 | sumOfRelaxedPrimalInfeasibilities_ = rhs.sumOfRelaxedPrimalInfeasibilities_; |
| 362 | infeasibilityWeight_ = rhs.infeasibilityWeight_; |
| 363 | numberDualInfeasibilities_ = rhs.numberDualInfeasibilities_; |
| 364 | numberPrimalInfeasibilities_ = rhs.numberPrimalInfeasibilities_; |
| 365 | noCheck_ = rhs.noCheck_; |
| 366 | firstGub_ = rhs.firstGub_; |
| 367 | lastGub_ = rhs.lastGub_; |
| 368 | gubType_ = rhs.gubType_; |
| 369 | model_ = rhs.model_; |
| 370 | } |
| 371 | return *this; |
412 | | // Assuming no gub rows deleted |
413 | | // We also assume all sets in same order |
414 | | // Get array with backward pointers |
415 | | int numberColumnsOld = rhs.matrix_->getNumCols(); |
416 | | int * array = new int [ numberColumnsOld]; |
417 | | int i; |
418 | | for (i = 0; i < numberColumnsOld; i++) |
419 | | array[i] = -1; |
420 | | for (int iSet = 0; iSet < numberSets_; iSet++) { |
421 | | for (int j = start_[iSet]; j < end_[iSet]; j++) |
422 | | array[j] = iSet; |
423 | | } |
424 | | numberSets_ = -1; |
425 | | int lastSet = -1; |
426 | | bool inSet = false; |
427 | | for (i = 0; i < numberColumns; i++) { |
428 | | int iColumn = whichColumns[i]; |
429 | | int iSet = array[iColumn]; |
430 | | if (iSet < 0) { |
431 | | inSet = false; |
432 | | } else { |
433 | | if (!inSet) { |
434 | | // start of new set but check okay |
435 | | if (iSet <= lastSet) |
436 | | throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix"); |
437 | | lastSet = iSet; |
438 | | numberSets_++; |
439 | | start_[numberSets_] = i; |
440 | | end_[numberSets_] = i + 1; |
441 | | lower_[numberSets_] = lower_[iSet]; |
442 | | upper_[numberSets_] = upper_[iSet]; |
443 | | inSet = true; |
444 | | } else { |
445 | | if (iSet < lastSet) { |
446 | | throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix"); |
447 | | } else if (iSet == lastSet) { |
448 | | end_[numberSets_] = i + 1; |
449 | | } else { |
450 | | // new set |
451 | | lastSet = iSet; |
452 | | numberSets_++; |
453 | | start_[numberSets_] = i; |
454 | | end_[numberSets_] = i + 1; |
455 | | lower_[numberSets_] = lower_[iSet]; |
456 | | upper_[numberSets_] = upper_[iSet]; |
457 | | } |
458 | | } |
459 | | } |
460 | | } |
461 | | delete [] array; |
462 | | numberSets_++; // adjust |
463 | | // Find type of gub |
464 | | firstGub_ = numberColumns + 1; |
465 | | lastGub_ = -1; |
466 | | for (i = 0; i < numberColumns; i++) { |
467 | | if (backward_[i] >= 0) { |
468 | | firstGub_ = CoinMin(firstGub_, i); |
469 | | lastGub_ = CoinMax(lastGub_, i); |
470 | | } |
471 | | } |
472 | | if (lastGub_ > 0) |
473 | | lastGub_++; |
474 | | gubType_ = 0; |
475 | | for (i = firstGub_; i < lastGub_; i++) { |
476 | | if (backward_[i] < 0) { |
477 | | gubType_ = 1; |
478 | | break; |
479 | | } |
480 | | } |
| 409 | // Assuming no gub rows deleted |
| 410 | // We also assume all sets in same order |
| 411 | // Get array with backward pointers |
| 412 | int numberColumnsOld = rhs.matrix_->getNumCols(); |
| 413 | int *array = new int[numberColumnsOld]; |
| 414 | int i; |
| 415 | for (i = 0; i < numberColumnsOld; i++) |
| 416 | array[i] = -1; |
| 417 | for (int iSet = 0; iSet < numberSets_; iSet++) { |
| 418 | for (int j = start_[iSet]; j < end_[iSet]; j++) |
| 419 | array[j] = iSet; |
| 420 | } |
| 421 | numberSets_ = -1; |
| 422 | int lastSet = -1; |
| 423 | bool inSet = false; |
| 424 | for (i = 0; i < numberColumns; i++) { |
| 425 | int iColumn = whichColumns[i]; |
| 426 | int iSet = array[iColumn]; |
| 427 | if (iSet < 0) { |
| 428 | inSet = false; |
| 429 | } else { |
| 430 | if (!inSet) { |
| 431 | // start of new set but check okay |
| 432 | if (iSet <= lastSet) |
| 433 | throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix"); |
| 434 | lastSet = iSet; |
| 435 | numberSets_++; |
| 436 | start_[numberSets_] = i; |
| 437 | end_[numberSets_] = i + 1; |
| 438 | lower_[numberSets_] = lower_[iSet]; |
| 439 | upper_[numberSets_] = upper_[iSet]; |
| 440 | inSet = true; |
| 441 | } else { |
| 442 | if (iSet < lastSet) { |
| 443 | throw CoinError("overlapping or non-monotonic sets", "subset constructor", "ClpGubMatrix"); |
| 444 | } else if (iSet == lastSet) { |
| 445 | end_[numberSets_] = i + 1; |
| 446 | } else { |
| 447 | // new set |
| 448 | lastSet = iSet; |
| 449 | numberSets_++; |
| 450 | start_[numberSets_] = i; |
| 451 | end_[numberSets_] = i + 1; |
| 452 | lower_[numberSets_] = lower_[iSet]; |
| 453 | upper_[numberSets_] = upper_[iSet]; |
| 454 | } |
| 455 | } |
| 456 | } |
| 457 | } |
| 458 | delete[] array; |
| 459 | numberSets_++; // adjust |
| 460 | // Find type of gub |
| 461 | firstGub_ = numberColumns + 1; |
| 462 | lastGub_ = -1; |
| 463 | for (i = 0; i < numberColumns; i++) { |
| 464 | if (backward_[i] >= 0) { |
| 465 | firstGub_ = CoinMin(firstGub_, i); |
| 466 | lastGub_ = CoinMax(lastGub_, i); |
| 467 | } |
| 468 | } |
| 469 | if (lastGub_ > 0) |
| 470 | lastGub_++; |
| 471 | gubType_ = 0; |
| 472 | for (i = firstGub_; i < lastGub_; i++) { |
| 473 | if (backward_[i] < 0) { |
| 474 | gubType_ = 1; |
| 475 | break; |
| 476 | } |
| 477 | } |
524 | | columnArray->clear(); |
525 | | double * pi = rowArray->denseVector(); |
526 | | int numberNonZero = 0; |
527 | | int * index = columnArray->getIndices(); |
528 | | double * array = columnArray->denseVector(); |
529 | | int numberInRowArray = rowArray->getNumElements(); |
530 | | // maybe I need one in OsiSimplex |
531 | | double zeroTolerance = model->zeroTolerance(); |
532 | | int numberRows = model->numberRows(); |
533 | | ClpPackedMatrix* rowCopy = |
534 | | dynamic_cast< ClpPackedMatrix*>(model->rowCopy()); |
535 | | bool packed = rowArray->packedMode(); |
536 | | double factor = 0.3; |
537 | | // We may not want to do by row if there may be cache problems |
538 | | int numberColumns = model->numberColumns(); |
539 | | // It would be nice to find L2 cache size - for moment 512K |
540 | | // Be slightly optimistic |
541 | | if (numberColumns * sizeof(double) > 1000000) { |
542 | | if (numberRows * 10 < numberColumns) |
543 | | factor = 0.1; |
544 | | else if (numberRows * 4 < numberColumns) |
545 | | factor = 0.15; |
546 | | else if (numberRows * 2 < numberColumns) |
547 | | factor = 0.2; |
548 | | //if (model->numberIterations()%50==0) |
549 | | //printf("%d nonzero\n",numberInRowArray); |
550 | | } |
551 | | // reduce for gub |
552 | | factor *= 0.5; |
553 | | assert (!y->getNumElements()); |
554 | | if (numberInRowArray > factor * numberRows || !rowCopy) { |
555 | | // do by column |
556 | | int iColumn; |
557 | | // get matrix data pointers |
558 | | const int * row = matrix_->getIndices(); |
559 | | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
560 | | const int * columnLength = matrix_->getVectorLengths(); |
561 | | const double * elementByColumn = matrix_->getElements(); |
562 | | const double * rowScale = model->rowScale(); |
563 | | int numberColumns = model->numberColumns(); |
564 | | int iSet = -1; |
565 | | double djMod = 0.0; |
566 | | if (packed) { |
567 | | // need to expand pi into y |
568 | | assert(y->capacity() >= numberRows); |
569 | | double * piOld = pi; |
570 | | pi = y->denseVector(); |
571 | | const int * whichRow = rowArray->getIndices(); |
572 | | int i; |
573 | | if (!rowScale) { |
574 | | // modify pi so can collapse to one loop |
575 | | for (i = 0; i < numberInRowArray; i++) { |
576 | | int iRow = whichRow[i]; |
577 | | pi[iRow] = scalar * piOld[i]; |
578 | | } |
579 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
580 | | if (backward_[iColumn] != iSet) { |
581 | | // get pi on gub row |
582 | | iSet = backward_[iColumn]; |
583 | | if (iSet >= 0) { |
584 | | int iBasic = keyVariable_[iSet]; |
585 | | if (iBasic < numberColumns) { |
586 | | // get dj without |
587 | | assert (model->getStatus(iBasic) == ClpSimplex::basic); |
588 | | djMod = 0.0; |
589 | | for (CoinBigIndex j = columnStart[iBasic]; |
590 | | j < columnStart[iBasic] + columnLength[iBasic]; j++) { |
591 | | int jRow = row[j]; |
592 | | djMod -= pi[jRow] * elementByColumn[j]; |
593 | | } |
594 | | } else { |
595 | | djMod = 0.0; |
596 | | } |
597 | | } else { |
598 | | djMod = 0.0; |
599 | | } |
600 | | } |
601 | | double value = -djMod; |
602 | | CoinBigIndex j; |
603 | | for (j = columnStart[iColumn]; |
604 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
605 | | int iRow = row[j]; |
606 | | value += pi[iRow] * elementByColumn[j]; |
607 | | } |
608 | | if (fabs(value) > zeroTolerance) { |
609 | | array[numberNonZero] = value; |
610 | | index[numberNonZero++] = iColumn; |
611 | | } |
612 | | } |
613 | | } else { |
614 | | // scaled |
615 | | // modify pi so can collapse to one loop |
616 | | for (i = 0; i < numberInRowArray; i++) { |
617 | | int iRow = whichRow[i]; |
618 | | pi[iRow] = scalar * piOld[i] * rowScale[iRow]; |
619 | | } |
620 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
621 | | if (backward_[iColumn] != iSet) { |
622 | | // get pi on gub row |
623 | | iSet = backward_[iColumn]; |
624 | | if (iSet >= 0) { |
625 | | int iBasic = keyVariable_[iSet]; |
626 | | if (iBasic < numberColumns) { |
627 | | // get dj without |
628 | | assert (model->getStatus(iBasic) == ClpSimplex::basic); |
629 | | djMod = 0.0; |
630 | | // scaled |
631 | | for (CoinBigIndex j = columnStart[iBasic]; |
632 | | j < columnStart[iBasic] + columnLength[iBasic]; j++) { |
633 | | int jRow = row[j]; |
634 | | djMod -= pi[jRow] * elementByColumn[j] * rowScale[jRow]; |
635 | | } |
636 | | } else { |
637 | | djMod = 0.0; |
638 | | } |
639 | | } else { |
640 | | djMod = 0.0; |
641 | | } |
642 | | } |
643 | | double value = -djMod; |
644 | | CoinBigIndex j; |
645 | | const double * columnScale = model->columnScale(); |
646 | | for (j = columnStart[iColumn]; |
647 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
648 | | int iRow = row[j]; |
649 | | value += pi[iRow] * elementByColumn[j]; |
650 | | } |
651 | | value *= columnScale[iColumn]; |
652 | | if (fabs(value) > zeroTolerance) { |
653 | | array[numberNonZero] = value; |
654 | | index[numberNonZero++] = iColumn; |
655 | | } |
656 | | } |
657 | | } |
658 | | // zero out |
659 | | for (i = 0; i < numberInRowArray; i++) { |
660 | | int iRow = whichRow[i]; |
661 | | pi[iRow] = 0.0; |
662 | | } |
663 | | } else { |
664 | | // code later |
665 | | assert (packed); |
666 | | if (!rowScale) { |
667 | | if (scalar == -1.0) { |
668 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
669 | | double value = 0.0; |
670 | | CoinBigIndex j; |
671 | | for (j = columnStart[iColumn]; |
672 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
673 | | int iRow = row[j]; |
674 | | value += pi[iRow] * elementByColumn[j]; |
675 | | } |
676 | | if (fabs(value) > zeroTolerance) { |
677 | | index[numberNonZero++] = iColumn; |
678 | | array[iColumn] = -value; |
679 | | } |
680 | | } |
681 | | } else if (scalar == 1.0) { |
682 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
683 | | double value = 0.0; |
684 | | CoinBigIndex j; |
685 | | for (j = columnStart[iColumn]; |
686 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
687 | | int iRow = row[j]; |
688 | | value += pi[iRow] * elementByColumn[j]; |
689 | | } |
690 | | if (fabs(value) > zeroTolerance) { |
691 | | index[numberNonZero++] = iColumn; |
692 | | array[iColumn] = value; |
693 | | } |
694 | | } |
695 | | } else { |
696 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
697 | | double value = 0.0; |
698 | | CoinBigIndex j; |
699 | | for (j = columnStart[iColumn]; |
700 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
701 | | int iRow = row[j]; |
702 | | value += pi[iRow] * elementByColumn[j]; |
703 | | } |
704 | | value *= scalar; |
705 | | if (fabs(value) > zeroTolerance) { |
706 | | index[numberNonZero++] = iColumn; |
707 | | array[iColumn] = value; |
708 | | } |
709 | | } |
710 | | } |
711 | | } else { |
712 | | // scaled |
713 | | if (scalar == -1.0) { |
714 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
715 | | double value = 0.0; |
716 | | CoinBigIndex j; |
717 | | const double * columnScale = model->columnScale(); |
718 | | for (j = columnStart[iColumn]; |
719 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
720 | | int iRow = row[j]; |
721 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
722 | | } |
723 | | value *= columnScale[iColumn]; |
724 | | if (fabs(value) > zeroTolerance) { |
725 | | index[numberNonZero++] = iColumn; |
726 | | array[iColumn] = -value; |
727 | | } |
728 | | } |
729 | | } else if (scalar == 1.0) { |
730 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
731 | | double value = 0.0; |
732 | | CoinBigIndex j; |
733 | | const double * columnScale = model->columnScale(); |
734 | | for (j = columnStart[iColumn]; |
735 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
736 | | int iRow = row[j]; |
737 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
738 | | } |
739 | | value *= columnScale[iColumn]; |
740 | | if (fabs(value) > zeroTolerance) { |
741 | | index[numberNonZero++] = iColumn; |
742 | | array[iColumn] = value; |
743 | | } |
744 | | } |
745 | | } else { |
746 | | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
747 | | double value = 0.0; |
748 | | CoinBigIndex j; |
749 | | const double * columnScale = model->columnScale(); |
750 | | for (j = columnStart[iColumn]; |
751 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
752 | | int iRow = row[j]; |
753 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
754 | | } |
755 | | value *= scalar * columnScale[iColumn]; |
756 | | if (fabs(value) > zeroTolerance) { |
757 | | index[numberNonZero++] = iColumn; |
758 | | array[iColumn] = value; |
759 | | } |
760 | | } |
761 | | } |
762 | | } |
763 | | } |
764 | | columnArray->setNumElements(numberNonZero); |
765 | | y->setNumElements(0); |
766 | | } else { |
767 | | // do by row |
768 | | transposeTimesByRow(model, scalar, rowArray, y, columnArray); |
769 | | } |
770 | | if (packed) |
771 | | columnArray->setPackedMode(true); |
772 | | if (0) { |
773 | | columnArray->checkClean(); |
774 | | int numberNonZero = columnArray->getNumElements();; |
775 | | int * index = columnArray->getIndices(); |
776 | | double * array = columnArray->denseVector(); |
777 | | int i; |
778 | | for (i = 0; i < numberNonZero; i++) { |
779 | | int j = index[i]; |
780 | | double value; |
781 | | if (packed) |
782 | | value = array[i]; |
783 | | else |
784 | | value = array[j]; |
785 | | printf("Ti %d %d %g\n", i, j, value); |
786 | | } |
787 | | } |
| 520 | columnArray->clear(); |
| 521 | double *pi = rowArray->denseVector(); |
| 522 | int numberNonZero = 0; |
| 523 | int *index = columnArray->getIndices(); |
| 524 | double *array = columnArray->denseVector(); |
| 525 | int numberInRowArray = rowArray->getNumElements(); |
| 526 | // maybe I need one in OsiSimplex |
| 527 | double zeroTolerance = model->zeroTolerance(); |
| 528 | int numberRows = model->numberRows(); |
| 529 | ClpPackedMatrix *rowCopy = dynamic_cast< ClpPackedMatrix * >(model->rowCopy()); |
| 530 | bool packed = rowArray->packedMode(); |
| 531 | double factor = 0.3; |
| 532 | // We may not want to do by row if there may be cache problems |
| 533 | int numberColumns = model->numberColumns(); |
| 534 | // It would be nice to find L2 cache size - for moment 512K |
| 535 | // Be slightly optimistic |
| 536 | if (numberColumns * sizeof(double) > 1000000) { |
| 537 | if (numberRows * 10 < numberColumns) |
| 538 | factor = 0.1; |
| 539 | else if (numberRows * 4 < numberColumns) |
| 540 | factor = 0.15; |
| 541 | else if (numberRows * 2 < numberColumns) |
| 542 | factor = 0.2; |
| 543 | //if (model->numberIterations()%50==0) |
| 544 | //printf("%d nonzero\n",numberInRowArray); |
| 545 | } |
| 546 | // reduce for gub |
| 547 | factor *= 0.5; |
| 548 | assert(!y->getNumElements()); |
| 549 | if (numberInRowArray > factor * numberRows || !rowCopy) { |
| 550 | // do by column |
| 551 | int iColumn; |
| 552 | // get matrix data pointers |
| 553 | const int *row = matrix_->getIndices(); |
| 554 | const CoinBigIndex *columnStart = matrix_->getVectorStarts(); |
| 555 | const int *columnLength = matrix_->getVectorLengths(); |
| 556 | const double *elementByColumn = matrix_->getElements(); |
| 557 | const double *rowScale = model->rowScale(); |
| 558 | int numberColumns = model->numberColumns(); |
| 559 | int iSet = -1; |
| 560 | double djMod = 0.0; |
| 561 | if (packed) { |
| 562 | // need to expand pi into y |
| 563 | assert(y->capacity() >= numberRows); |
| 564 | double *piOld = pi; |
| 565 | pi = y->denseVector(); |
| 566 | const int *whichRow = rowArray->getIndices(); |
| 567 | int i; |
| 568 | if (!rowScale) { |
| 569 | // modify pi so can collapse to one loop |
| 570 | for (i = 0; i < numberInRowArray; i++) { |
| 571 | int iRow = whichRow[i]; |
| 572 | pi[iRow] = scalar * piOld[i]; |
| 573 | } |
| 574 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 575 | if (backward_[iColumn] != iSet) { |
| 576 | // get pi on gub row |
| 577 | iSet = backward_[iColumn]; |
| 578 | if (iSet >= 0) { |
| 579 | int iBasic = keyVariable_[iSet]; |
| 580 | if (iBasic < numberColumns) { |
| 581 | // get dj without |
| 582 | assert(model->getStatus(iBasic) == ClpSimplex::basic); |
| 583 | djMod = 0.0; |
| 584 | for (CoinBigIndex j = columnStart[iBasic]; |
| 585 | j < columnStart[iBasic] + columnLength[iBasic]; j++) { |
| 586 | int jRow = row[j]; |
| 587 | djMod -= pi[jRow] * elementByColumn[j]; |
| 588 | } |
| 589 | } else { |
| 590 | djMod = 0.0; |
| 591 | } |
| 592 | } else { |
| 593 | djMod = 0.0; |
| 594 | } |
| 595 | } |
| 596 | double value = -djMod; |
| 597 | CoinBigIndex j; |
| 598 | for (j = columnStart[iColumn]; |
| 599 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 600 | int iRow = row[j]; |
| 601 | value += pi[iRow] * elementByColumn[j]; |
| 602 | } |
| 603 | if (fabs(value) > zeroTolerance) { |
| 604 | array[numberNonZero] = value; |
| 605 | index[numberNonZero++] = iColumn; |
| 606 | } |
| 607 | } |
| 608 | } else { |
| 609 | // scaled |
| 610 | // modify pi so can collapse to one loop |
| 611 | for (i = 0; i < numberInRowArray; i++) { |
| 612 | int iRow = whichRow[i]; |
| 613 | pi[iRow] = scalar * piOld[i] * rowScale[iRow]; |
| 614 | } |
| 615 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 616 | if (backward_[iColumn] != iSet) { |
| 617 | // get pi on gub row |
| 618 | iSet = backward_[iColumn]; |
| 619 | if (iSet >= 0) { |
| 620 | int iBasic = keyVariable_[iSet]; |
| 621 | if (iBasic < numberColumns) { |
| 622 | // get dj without |
| 623 | assert(model->getStatus(iBasic) == ClpSimplex::basic); |
| 624 | djMod = 0.0; |
| 625 | // scaled |
| 626 | for (CoinBigIndex j = columnStart[iBasic]; |
| 627 | j < columnStart[iBasic] + columnLength[iBasic]; j++) { |
| 628 | int jRow = row[j]; |
| 629 | djMod -= pi[jRow] * elementByColumn[j] * rowScale[jRow]; |
| 630 | } |
| 631 | } else { |
| 632 | djMod = 0.0; |
| 633 | } |
| 634 | } else { |
| 635 | djMod = 0.0; |
| 636 | } |
| 637 | } |
| 638 | double value = -djMod; |
| 639 | CoinBigIndex j; |
| 640 | const double *columnScale = model->columnScale(); |
| 641 | for (j = columnStart[iColumn]; |
| 642 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 643 | int iRow = row[j]; |
| 644 | value += pi[iRow] * elementByColumn[j]; |
| 645 | } |
| 646 | value *= columnScale[iColumn]; |
| 647 | if (fabs(value) > zeroTolerance) { |
| 648 | array[numberNonZero] = value; |
| 649 | index[numberNonZero++] = iColumn; |
| 650 | } |
| 651 | } |
| 652 | } |
| 653 | // zero out |
| 654 | for (i = 0; i < numberInRowArray; i++) { |
| 655 | int iRow = whichRow[i]; |
| 656 | pi[iRow] = 0.0; |
| 657 | } |
| 658 | } else { |
| 659 | // code later |
| 660 | assert(packed); |
| 661 | if (!rowScale) { |
| 662 | if (scalar == -1.0) { |
| 663 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 664 | double value = 0.0; |
| 665 | CoinBigIndex j; |
| 666 | for (j = columnStart[iColumn]; |
| 667 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 668 | int iRow = row[j]; |
| 669 | value += pi[iRow] * elementByColumn[j]; |
| 670 | } |
| 671 | if (fabs(value) > zeroTolerance) { |
| 672 | index[numberNonZero++] = iColumn; |
| 673 | array[iColumn] = -value; |
| 674 | } |
| 675 | } |
| 676 | } else if (scalar == 1.0) { |
| 677 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 678 | double value = 0.0; |
| 679 | CoinBigIndex j; |
| 680 | for (j = columnStart[iColumn]; |
| 681 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 682 | int iRow = row[j]; |
| 683 | value += pi[iRow] * elementByColumn[j]; |
| 684 | } |
| 685 | if (fabs(value) > zeroTolerance) { |
| 686 | index[numberNonZero++] = iColumn; |
| 687 | array[iColumn] = value; |
| 688 | } |
| 689 | } |
| 690 | } else { |
| 691 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 692 | double value = 0.0; |
| 693 | CoinBigIndex j; |
| 694 | for (j = columnStart[iColumn]; |
| 695 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 696 | int iRow = row[j]; |
| 697 | value += pi[iRow] * elementByColumn[j]; |
| 698 | } |
| 699 | value *= scalar; |
| 700 | if (fabs(value) > zeroTolerance) { |
| 701 | index[numberNonZero++] = iColumn; |
| 702 | array[iColumn] = value; |
| 703 | } |
| 704 | } |
| 705 | } |
| 706 | } else { |
| 707 | // scaled |
| 708 | if (scalar == -1.0) { |
| 709 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 710 | double value = 0.0; |
| 711 | CoinBigIndex j; |
| 712 | const double *columnScale = model->columnScale(); |
| 713 | for (j = columnStart[iColumn]; |
| 714 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 715 | int iRow = row[j]; |
| 716 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 717 | } |
| 718 | value *= columnScale[iColumn]; |
| 719 | if (fabs(value) > zeroTolerance) { |
| 720 | index[numberNonZero++] = iColumn; |
| 721 | array[iColumn] = -value; |
| 722 | } |
| 723 | } |
| 724 | } else if (scalar == 1.0) { |
| 725 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 726 | double value = 0.0; |
| 727 | CoinBigIndex j; |
| 728 | const double *columnScale = model->columnScale(); |
| 729 | for (j = columnStart[iColumn]; |
| 730 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 731 | int iRow = row[j]; |
| 732 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 733 | } |
| 734 | value *= columnScale[iColumn]; |
| 735 | if (fabs(value) > zeroTolerance) { |
| 736 | index[numberNonZero++] = iColumn; |
| 737 | array[iColumn] = value; |
| 738 | } |
| 739 | } |
| 740 | } else { |
| 741 | for (iColumn = 0; iColumn < numberColumns; iColumn++) { |
| 742 | double value = 0.0; |
| 743 | CoinBigIndex j; |
| 744 | const double *columnScale = model->columnScale(); |
| 745 | for (j = columnStart[iColumn]; |
| 746 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 747 | int iRow = row[j]; |
| 748 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 749 | } |
| 750 | value *= scalar * columnScale[iColumn]; |
| 751 | if (fabs(value) > zeroTolerance) { |
| 752 | index[numberNonZero++] = iColumn; |
| 753 | array[iColumn] = value; |
| 754 | } |
| 755 | } |
| 756 | } |
| 757 | } |
| 758 | } |
| 759 | columnArray->setNumElements(numberNonZero); |
| 760 | y->setNumElements(0); |
| 761 | } else { |
| 762 | // do by row |
| 763 | transposeTimesByRow(model, scalar, rowArray, y, columnArray); |
| 764 | } |
| 765 | if (packed) |
| 766 | columnArray->setPackedMode(true); |
| 767 | if (0) { |
| 768 | columnArray->checkClean(); |
| 769 | int numberNonZero = columnArray->getNumElements(); |
| 770 | ; |
| 771 | int *index = columnArray->getIndices(); |
| 772 | double *array = columnArray->denseVector(); |
| 773 | int i; |
| 774 | for (i = 0; i < numberNonZero; i++) { |
| 775 | int j = index[i]; |
| 776 | double value; |
| 777 | if (packed) |
| 778 | value = array[i]; |
| 779 | else |
| 780 | value = array[j]; |
| 781 | printf("Ti %d %d %g\n", i, j, value); |
| 782 | } |
| 783 | } |
813 | | columnArray->clear(); |
814 | | double * pi = rowArray->denseVector(); |
815 | | double * array = columnArray->denseVector(); |
816 | | int jColumn; |
817 | | // get matrix data pointers |
818 | | const int * row = matrix_->getIndices(); |
819 | | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
820 | | const int * columnLength = matrix_->getVectorLengths(); |
821 | | const double * elementByColumn = matrix_->getElements(); |
822 | | const double * rowScale = model->rowScale(); |
823 | | int numberToDo = y->getNumElements(); |
824 | | const int * which = y->getIndices(); |
825 | | assert (!rowArray->packedMode()); |
826 | | columnArray->setPacked(); |
827 | | int numberTouched = 0; |
828 | | if (!rowScale) { |
829 | | for (jColumn = 0; jColumn < numberToDo; jColumn++) { |
830 | | int iColumn = which[jColumn]; |
831 | | double value = 0.0; |
832 | | CoinBigIndex j; |
833 | | for (j = columnStart[iColumn]; |
834 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
835 | | int iRow = row[j]; |
836 | | value += pi[iRow] * elementByColumn[j]; |
837 | | } |
838 | | array[jColumn] = value; |
839 | | if (value) { |
840 | | int iSet = backward_[iColumn]; |
841 | | if (iSet >= 0) { |
842 | | int iBasic = keyVariable_[iSet]; |
843 | | if (iBasic == iColumn) { |
844 | | toIndex_[iSet] = jColumn; |
845 | | fromIndex_[numberTouched++] = iSet; |
846 | | } |
847 | | } |
848 | | } |
849 | | } |
850 | | } else { |
851 | | // scaled |
852 | | for (jColumn = 0; jColumn < numberToDo; jColumn++) { |
853 | | int iColumn = which[jColumn]; |
854 | | double value = 0.0; |
855 | | CoinBigIndex j; |
856 | | const double * columnScale = model->columnScale(); |
857 | | for (j = columnStart[iColumn]; |
858 | | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
859 | | int iRow = row[j]; |
860 | | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
861 | | } |
862 | | value *= columnScale[iColumn]; |
863 | | array[jColumn] = value; |
864 | | if (value) { |
865 | | int iSet = backward_[iColumn]; |
866 | | if (iSet >= 0) { |
867 | | int iBasic = keyVariable_[iSet]; |
868 | | if (iBasic == iColumn) { |
869 | | toIndex_[iSet] = jColumn; |
870 | | fromIndex_[numberTouched++] = iSet; |
871 | | } |
872 | | } |
873 | | } |
874 | | } |
875 | | } |
876 | | // adjust djs |
877 | | for (jColumn = 0; jColumn < numberToDo; jColumn++) { |
878 | | int iColumn = which[jColumn]; |
879 | | int iSet = backward_[iColumn]; |
880 | | if (iSet >= 0) { |
881 | | int kColumn = toIndex_[iSet]; |
882 | | if (kColumn >= 0) |
883 | | array[jColumn] -= array[kColumn]; |
884 | | } |
885 | | } |
886 | | // and clear basic |
887 | | for (int j = 0; j < numberTouched; j++) { |
888 | | int iSet = fromIndex_[j]; |
889 | | int kColumn = toIndex_[iSet]; |
890 | | toIndex_[iSet] = -1; |
891 | | array[kColumn] = 0.0; |
892 | | } |
| 807 | columnArray->clear(); |
| 808 | double *pi = rowArray->denseVector(); |
| 809 | double *array = columnArray->denseVector(); |
| 810 | int jColumn; |
| 811 | // get matrix data pointers |
| 812 | const int *row = matrix_->getIndices(); |
| 813 | const CoinBigIndex *columnStart = matrix_->getVectorStarts(); |
| 814 | const int *columnLength = matrix_->getVectorLengths(); |
| 815 | const double *elementByColumn = matrix_->getElements(); |
| 816 | const double *rowScale = model->rowScale(); |
| 817 | int numberToDo = y->getNumElements(); |
| 818 | const int *which = y->getIndices(); |
| 819 | assert(!rowArray->packedMode()); |
| 820 | columnArray->setPacked(); |
| 821 | int numberTouched = 0; |
| 822 | if (!rowScale) { |
| 823 | for (jColumn = 0; jColumn < numberToDo; jColumn++) { |
| 824 | int iColumn = which[jColumn]; |
| 825 | double value = 0.0; |
| 826 | CoinBigIndex j; |
| 827 | for (j = columnStart[iColumn]; |
| 828 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 829 | int iRow = row[j]; |
| 830 | value += pi[iRow] * elementByColumn[j]; |
| 831 | } |
| 832 | array[jColumn] = value; |
| 833 | if (value) { |
| 834 | int iSet = backward_[iColumn]; |
| 835 | if (iSet >= 0) { |
| 836 | int iBasic = keyVariable_[iSet]; |
| 837 | if (iBasic == iColumn) { |
| 838 | toIndex_[iSet] = jColumn; |
| 839 | fromIndex_[numberTouched++] = iSet; |
| 840 | } |
| 841 | } |
| 842 | } |
| 843 | } |
| 844 | } else { |
| 845 | // scaled |
| 846 | for (jColumn = 0; jColumn < numberToDo; jColumn++) { |
| 847 | int iColumn = which[jColumn]; |
| 848 | double value = 0.0; |
| 849 | CoinBigIndex j; |
| 850 | const double *columnScale = model->columnScale(); |
| 851 | for (j = columnStart[iColumn]; |
| 852 | j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 853 | int iRow = row[j]; |
| 854 | value += pi[iRow] * elementByColumn[j] * rowScale[iRow]; |
| 855 | } |
| 856 | value *= columnScale[iColumn]; |
| 857 | array[jColumn] = value; |
| 858 | if (value) { |
| 859 | int iSet = backward_[iColumn]; |
| 860 | if (iSet >= 0) { |
| 861 | int iBasic = keyVariable_[iSet]; |
| 862 | if (iBasic == iColumn) { |
| 863 | toIndex_[iSet] = jColumn; |
| 864 | fromIndex_[numberTouched++] = iSet; |
| 865 | } |
| 866 | } |
| 867 | } |
| 868 | } |
| 869 | } |
| 870 | // adjust djs |
| 871 | for (jColumn = 0; jColumn < numberToDo; jColumn++) { |
| 872 | int iColumn = which[jColumn]; |
| 873 | int iSet = backward_[iColumn]; |
| 874 | if (iSet >= 0) { |
| 875 | int kColumn = toIndex_[iSet]; |
| 876 | if (kColumn >= 0) |
| 877 | array[jColumn] -= array[kColumn]; |
| 878 | } |
| 879 | } |
| 880 | // and clear basic |
| 881 | for (int j = 0; j < numberTouched; j++) { |
| 882 | int iSet = fromIndex_[j]; |
| 883 | int kColumn = toIndex_[iSet]; |
| 884 | toIndex_[iSet] = -1; |
| 885 | array[kColumn] = 0.0; |
| 886 | } |
984 | | int i; |
985 | | int numberColumns = getNumCols(); |
986 | | const int * columnLength = matrix_->getVectorLengths(); |
987 | | int numberRows = getNumRows(); |
988 | | assert (next_ || !elementU) ; |
989 | | CoinBigIndex numberElements = start[0]; |
990 | | int lastSet = -1; |
991 | | int key = -1; |
992 | | int keyLength = -1; |
993 | | double * work = new double[numberRows]; |
994 | | CoinZeroN(work, numberRows); |
995 | | char * mark = new char[numberRows]; |
996 | | CoinZeroN(mark, numberRows); |
997 | | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
998 | | const int * row = matrix_->getIndices(); |
999 | | const double * elementByColumn = matrix_->getElements(); |
1000 | | const double * rowScale = model->rowScale(); |
1001 | | int numberBasic = 0; |
1002 | | if (0) { |
1003 | | printf("%d basiccolumns\n", numberColumnBasic); |
1004 | | int i; |
1005 | | for (i = 0; i < numberSets_; i++) { |
1006 | | int k = keyVariable_[i]; |
1007 | | if (k < numberColumns) { |
1008 | | printf("key %d on set %d, %d elements\n", k, i, columnStart[k+1] - columnStart[k]); |
1009 | | for (int j = columnStart[k]; j < columnStart[k+1]; j++) |
1010 | | printf("row %d el %g\n", row[j], elementByColumn[j]); |
1011 | | } else { |
1012 | | printf("slack key on set %d\n", i); |
1013 | | } |
1014 | | } |
1015 | | } |
1016 | | // fill |
1017 | | if (!rowScale) { |
1018 | | // no scaling |
1019 | | for (i = 0; i < numberColumnBasic; i++) { |
1020 | | int iColumn = whichColumn[i]; |
1021 | | int iSet = backward_[iColumn]; |
1022 | | int length = columnLength[iColumn]; |
1023 | | if (0) { |
1024 | | int k = iColumn; |
1025 | | printf("column %d in set %d, %d elements\n", k, iSet, columnStart[k+1] - columnStart[k]); |
1026 | | for (int j = columnStart[k]; j < columnStart[k+1]; j++) |
1027 | | printf("row %d el %g\n", row[j], elementByColumn[j]); |
1028 | | } |
1029 | | CoinBigIndex j; |
1030 | | if (iSet < 0 || keyVariable_[iSet] >= numberColumns) { |
1031 | | for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
1032 | | double value = elementByColumn[j]; |
1033 | | if (fabs(value) > 1.0e-20) { |
1034 | | int iRow = row[j]; |
1035 | | indexRowU[numberElements] = iRow; |
1036 | | rowCount[iRow]++; |
1037 | | elementU[numberElements++] = value; |
1038 | | } |
1039 | | } |
1040 | | // end of column |
1041 | | columnCount[numberBasic] = numberElements - start[numberBasic]; |
1042 | | numberBasic++; |
1043 | | start[numberBasic] = numberElements; |
1044 | | } else { |
1045 | | // in gub set |
1046 | | if (iColumn != keyVariable_[iSet]) { |
1047 | | // not key |
1048 | | if (lastSet != iSet) { |
1049 | | // erase work |
1050 | | if (key >= 0) { |
1051 | | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
1052 | | int iRow = row[j]; |
1053 | | work[iRow] = 0.0; |
1054 | | mark[iRow] = 0; |
1055 | | } |
1056 | | } |
1057 | | key = keyVariable_[iSet]; |
1058 | | lastSet = iSet; |
1059 | | keyLength = columnLength[key]; |
1060 | | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
1061 | | int iRow = row[j]; |
1062 | | work[iRow] = elementByColumn[j]; |
1063 | | mark[iRow] = 1; |
1064 | | } |
1065 | | } |
1066 | | for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) { |
1067 | | int iRow = row[j]; |
1068 | | double value = elementByColumn[j]; |
1069 | | if (mark[iRow]) { |
1070 | | mark[iRow] = 0; |
1071 | | double keyValue = work[iRow]; |
1072 | | value -= keyValue; |
1073 | | } |
1074 | | if (fabs(value) > 1.0e-20) { |
1075 | | indexRowU[numberElements] = iRow; |
1076 | | rowCount[iRow]++; |
1077 | | elementU[numberElements++] = value; |
1078 | | } |
1079 | | } |
1080 | | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
1081 | | int iRow = row[j]; |
1082 | | if (mark[iRow]) { |
1083 | | double value = -work[iRow]; |
1084 | | if (fabs(value) > 1.0e-20) { |
1085 | | indexRowU[numberElements] = iRow; |
1086 | | rowCount[iRow]++; |
1087 | | elementU[numberElements++] = value; |
1088 | | } |
1089 | | } else { |
1090 | | // just put back mark |
1091 | | mark[iRow] = 1; |
1092 | | } |
1093 | | } |
1094 | | // end of column |
1095 | | columnCount[numberBasic] = numberElements - start[numberBasic]; |
1096 | | numberBasic++; |
1097 | | start[numberBasic] = numberElements; |
1098 | | } |
1099 | | } |
1100 | | } |
1101 | | } else { |
1102 | | // scaling |
1103 | | const double * columnScale = model->columnScale(); |
1104 | | for (i = 0; i < numberColumnBasic; i++) { |
1105 | | int iColumn = whichColumn[i]; |
1106 | | int iSet = backward_[iColumn]; |
1107 | | int length = columnLength[iColumn]; |
1108 | | CoinBigIndex j; |
1109 | | if (iSet < 0 || keyVariable_[iSet] >= numberColumns) { |
1110 | | double scale = columnScale[iColumn]; |
1111 | | for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
1112 | | int iRow = row[j]; |
1113 | | double value = elementByColumn[j] * scale * rowScale[iRow]; |
1114 | | if (fabs(value) > 1.0e-20) { |
1115 | | indexRowU[numberElements] = iRow; |
1116 | | rowCount[iRow]++; |
1117 | | elementU[numberElements++] = value; |
1118 | | } |
1119 | | } |
1120 | | // end of column |
1121 | | columnCount[numberBasic] = numberElements - start[numberBasic]; |
1122 | | numberBasic++; |
1123 | | start[numberBasic] = numberElements; |
1124 | | } else { |
1125 | | // in gub set |
1126 | | if (iColumn != keyVariable_[iSet]) { |
1127 | | double scale = columnScale[iColumn]; |
1128 | | // not key |
1129 | | if (lastSet < iSet) { |
1130 | | // erase work |
1131 | | if (key >= 0) { |
1132 | | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
1133 | | int iRow = row[j]; |
1134 | | work[iRow] = 0.0; |
1135 | | mark[iRow] = 0; |
1136 | | } |
1137 | | } |
1138 | | key = keyVariable_[iSet]; |
1139 | | lastSet = iSet; |
1140 | | keyLength = columnLength[key]; |
1141 | | double scale = columnScale[key]; |
1142 | | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
1143 | | int iRow = row[j]; |
1144 | | work[iRow] = elementByColumn[j] * scale * rowScale[iRow]; |
1145 | | mark[iRow] = 1; |
1146 | | } |
1147 | | } |
1148 | | for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) { |
1149 | | int iRow = row[j]; |
1150 | | double value = elementByColumn[j] * scale * rowScale[iRow]; |
1151 | | if (mark[iRow]) { |
1152 | | mark[iRow] = 0; |
1153 | | double keyValue = work[iRow]; |
1154 | | value -= keyValue; |
1155 | | } |
1156 | | if (fabs(value) > 1.0e-20) { |
1157 | | indexRowU[numberElements] = iRow; |
1158 | | rowCount[iRow]++; |
1159 | | elementU[numberElements++] = value; |
1160 | | } |
1161 | | } |
1162 | | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
1163 | | int iRow = row[j]; |
1164 | | if (mark[iRow]) { |
1165 | | double value = -work[iRow]; |
1166 | | if (fabs(value) > 1.0e-20) { |
1167 | | indexRowU[numberElements] = iRow; |
1168 | | rowCount[iRow]++; |
1169 | | elementU[numberElements++] = value; |
1170 | | } |
1171 | | } else { |
1172 | | // just put back mark |
1173 | | mark[iRow] = 1; |
1174 | | } |
1175 | | } |
1176 | | // end of column |
1177 | | columnCount[numberBasic] = numberElements - start[numberBasic]; |
1178 | | numberBasic++; |
1179 | | start[numberBasic] = numberElements; |
1180 | | } |
1181 | | } |
1182 | | } |
1183 | | } |
1184 | | delete [] work; |
1185 | | delete [] mark; |
1186 | | // update number of column basic |
1187 | | numberColumnBasic = numberBasic; |
| 976 | int i; |
| 977 | int numberColumns = getNumCols(); |
| 978 | const int *columnLength = matrix_->getVectorLengths(); |
| 979 | int numberRows = getNumRows(); |
| 980 | assert(next_ || !elementU); |
| 981 | CoinBigIndex numberElements = start[0]; |
| 982 | int lastSet = -1; |
| 983 | int key = -1; |
| 984 | int keyLength = -1; |
| 985 | double *work = new double[numberRows]; |
| 986 | CoinZeroN(work, numberRows); |
| 987 | char *mark = new char[numberRows]; |
| 988 | CoinZeroN(mark, numberRows); |
| 989 | const CoinBigIndex *columnStart = matrix_->getVectorStarts(); |
| 990 | const int *row = matrix_->getIndices(); |
| 991 | const double *elementByColumn = matrix_->getElements(); |
| 992 | const double *rowScale = model->rowScale(); |
| 993 | int numberBasic = 0; |
| 994 | if (0) { |
| 995 | printf("%d basiccolumns\n", numberColumnBasic); |
| 996 | int i; |
| 997 | for (i = 0; i < numberSets_; i++) { |
| 998 | int k = keyVariable_[i]; |
| 999 | if (k < numberColumns) { |
| 1000 | printf("key %d on set %d, %d elements\n", k, i, columnStart[k + 1] - columnStart[k]); |
| 1001 | for (int j = columnStart[k]; j < columnStart[k + 1]; j++) |
| 1002 | printf("row %d el %g\n", row[j], elementByColumn[j]); |
| 1003 | } else { |
| 1004 | printf("slack key on set %d\n", i); |
| 1005 | } |
| 1006 | } |
| 1007 | } |
| 1008 | // fill |
| 1009 | if (!rowScale) { |
| 1010 | // no scaling |
| 1011 | for (i = 0; i < numberColumnBasic; i++) { |
| 1012 | int iColumn = whichColumn[i]; |
| 1013 | int iSet = backward_[iColumn]; |
| 1014 | int length = columnLength[iColumn]; |
| 1015 | if (0) { |
| 1016 | int k = iColumn; |
| 1017 | printf("column %d in set %d, %d elements\n", k, iSet, columnStart[k + 1] - columnStart[k]); |
| 1018 | for (int j = columnStart[k]; j < columnStart[k + 1]; j++) |
| 1019 | printf("row %d el %g\n", row[j], elementByColumn[j]); |
| 1020 | } |
| 1021 | CoinBigIndex j; |
| 1022 | if (iSet < 0 || keyVariable_[iSet] >= numberColumns) { |
| 1023 | for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 1024 | double value = elementByColumn[j]; |
| 1025 | if (fabs(value) > 1.0e-20) { |
| 1026 | int iRow = row[j]; |
| 1027 | indexRowU[numberElements] = iRow; |
| 1028 | rowCount[iRow]++; |
| 1029 | elementU[numberElements++] = value; |
| 1030 | } |
| 1031 | } |
| 1032 | // end of column |
| 1033 | columnCount[numberBasic] = numberElements - start[numberBasic]; |
| 1034 | numberBasic++; |
| 1035 | start[numberBasic] = numberElements; |
| 1036 | } else { |
| 1037 | // in gub set |
| 1038 | if (iColumn != keyVariable_[iSet]) { |
| 1039 | // not key |
| 1040 | if (lastSet != iSet) { |
| 1041 | // erase work |
| 1042 | if (key >= 0) { |
| 1043 | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
| 1044 | int iRow = row[j]; |
| 1045 | work[iRow] = 0.0; |
| 1046 | mark[iRow] = 0; |
| 1047 | } |
| 1048 | } |
| 1049 | key = keyVariable_[iSet]; |
| 1050 | lastSet = iSet; |
| 1051 | keyLength = columnLength[key]; |
| 1052 | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
| 1053 | int iRow = row[j]; |
| 1054 | work[iRow] = elementByColumn[j]; |
| 1055 | mark[iRow] = 1; |
| 1056 | } |
| 1057 | } |
| 1058 | for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) { |
| 1059 | int iRow = row[j]; |
| 1060 | double value = elementByColumn[j]; |
| 1061 | if (mark[iRow]) { |
| 1062 | mark[iRow] = 0; |
| 1063 | double keyValue = work[iRow]; |
| 1064 | value -= keyValue; |
| 1065 | } |
| 1066 | if (fabs(value) > 1.0e-20) { |
| 1067 | indexRowU[numberElements] = iRow; |
| 1068 | rowCount[iRow]++; |
| 1069 | elementU[numberElements++] = value; |
| 1070 | } |
| 1071 | } |
| 1072 | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
| 1073 | int iRow = row[j]; |
| 1074 | if (mark[iRow]) { |
| 1075 | double value = -work[iRow]; |
| 1076 | if (fabs(value) > 1.0e-20) { |
| 1077 | indexRowU[numberElements] = iRow; |
| 1078 | rowCount[iRow]++; |
| 1079 | elementU[numberElements++] = value; |
| 1080 | } |
| 1081 | } else { |
| 1082 | // just put back mark |
| 1083 | mark[iRow] = 1; |
| 1084 | } |
| 1085 | } |
| 1086 | // end of column |
| 1087 | columnCount[numberBasic] = numberElements - start[numberBasic]; |
| 1088 | numberBasic++; |
| 1089 | start[numberBasic] = numberElements; |
| 1090 | } |
| 1091 | } |
| 1092 | } |
| 1093 | } else { |
| 1094 | // scaling |
| 1095 | const double *columnScale = model->columnScale(); |
| 1096 | for (i = 0; i < numberColumnBasic; i++) { |
| 1097 | int iColumn = whichColumn[i]; |
| 1098 | int iSet = backward_[iColumn]; |
| 1099 | int length = columnLength[iColumn]; |
| 1100 | CoinBigIndex j; |
| 1101 | if (iSet < 0 || keyVariable_[iSet] >= numberColumns) { |
| 1102 | double scale = columnScale[iColumn]; |
| 1103 | for (j = columnStart[iColumn]; j < columnStart[iColumn] + columnLength[iColumn]; j++) { |
| 1104 | int iRow = row[j]; |
| 1105 | double value = elementByColumn[j] * scale * rowScale[iRow]; |
| 1106 | if (fabs(value) > 1.0e-20) { |
| 1107 | indexRowU[numberElements] = iRow; |
| 1108 | rowCount[iRow]++; |
| 1109 | elementU[numberElements++] = value; |
| 1110 | } |
| 1111 | } |
| 1112 | // end of column |
| 1113 | columnCount[numberBasic] = numberElements - start[numberBasic]; |
| 1114 | numberBasic++; |
| 1115 | start[numberBasic] = numberElements; |
| 1116 | } else { |
| 1117 | // in gub set |
| 1118 | if (iColumn != keyVariable_[iSet]) { |
| 1119 | double scale = columnScale[iColumn]; |
| 1120 | // not key |
| 1121 | if (lastSet < iSet) { |
| 1122 | // erase work |
| 1123 | if (key >= 0) { |
| 1124 | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
| 1125 | int iRow = row[j]; |
| 1126 | work[iRow] = 0.0; |
| 1127 | mark[iRow] = 0; |
| 1128 | } |
| 1129 | } |
| 1130 | key = keyVariable_[iSet]; |
| 1131 | lastSet = iSet; |
| 1132 | keyLength = columnLength[key]; |
| 1133 | double scale = columnScale[key]; |
| 1134 | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
| 1135 | int iRow = row[j]; |
| 1136 | work[iRow] = elementByColumn[j] * scale * rowScale[iRow]; |
| 1137 | mark[iRow] = 1; |
| 1138 | } |
| 1139 | } |
| 1140 | for (j = columnStart[iColumn]; j < columnStart[iColumn] + length; j++) { |
| 1141 | int iRow = row[j]; |
| 1142 | double value = elementByColumn[j] * scale * rowScale[iRow]; |
| 1143 | if (mark[iRow]) { |
| 1144 | mark[iRow] = 0; |
| 1145 | double keyValue = work[iRow]; |
| 1146 | value -= keyValue; |
| 1147 | } |
| 1148 | if (fabs(value) > 1.0e-20) { |
| 1149 | indexRowU[numberElements] = iRow; |
| 1150 | rowCount[iRow]++; |
| 1151 | elementU[numberElements++] = value; |
| 1152 | } |
| 1153 | } |
| 1154 | for (j = columnStart[key]; j < columnStart[key] + keyLength; j++) { |
| 1155 | int iRow = row[j]; |
| 1156 | if (mark[iRow]) { |
| 1157 | double value = -work[iRow]; |
| 1158 | if (fabs(value) > 1.0e-20) { |
| 1159 | indexRowU[numberElements] = iRow; |
| 1160 | rowCount[iRow]++; |
| 1161 | elementU[numberElements++] = value; |
| 1162 | } |
| 1163 | } else { |
| 1164 | // just put back mark |
| 1165 | mark[iRow] = 1; |
| 1166 | } |
| 1167 | } |
| 1168 | // end of column |
| 1169 | columnCount[numberBasic] = numberElements - start[numberBasic]; |
| 1170 | numberBasic++; |
| 1171 | start[numberBasic] = numberElements; |
| 1172 | } |
| 1173 | } |
| 1174 | } |
| 1175 | } |
| 1176 | delete[] work; |
| 1177 | delete[] mark; |
| 1178 | // update number of column basic |
| 1179 | numberColumnBasic = numberBasic; |
1215 | | int numberColumns = model->numberColumns(); |
1216 | | if (iColumn < numberColumns) { |
1217 | | // Do packed part |
1218 | | ClpPackedMatrix::unpackPacked(model, rowArray, iColumn); |
1219 | | int iSet = backward_[iColumn]; |
1220 | | if (iSet >= 0) { |
1221 | | // columns are in order |
1222 | | int iBasic = keyVariable_[iSet]; |
1223 | | if (iBasic < numberColumns) { |
1224 | | int number = rowArray->getNumElements(); |
1225 | | const double * rowScale = model->rowScale(); |
1226 | | const int * row = matrix_->getIndices(); |
1227 | | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
1228 | | const int * columnLength = matrix_->getVectorLengths(); |
1229 | | const double * elementByColumn = matrix_->getElements(); |
1230 | | double * array = rowArray->denseVector(); |
1231 | | int * index = rowArray->getIndices(); |
1232 | | CoinBigIndex i; |
1233 | | int numberOld = number; |
1234 | | int lastIndex = 0; |
1235 | | int next = index[lastIndex]; |
1236 | | if (!rowScale) { |
1237 | | for (i = columnStart[iBasic]; |
1238 | | i < columnStart[iBasic] + columnLength[iBasic]; i++) { |
1239 | | int iRow = row[i]; |
1240 | | while (iRow > next) { |
1241 | | lastIndex++; |
1242 | | if (lastIndex == numberOld) |
1243 | | next = matrix_->getNumRows(); |
1244 | | else |
1245 | | next = index[lastIndex]; |
1246 | | } |
1247 | | if (iRow < next) { |
1248 | | array[number] = -elementByColumn[i]; |
1249 | | index[number++] = iRow; |
1250 | | } else { |
1251 | | assert (iRow == next); |
1252 | | array[lastIndex] -= elementByColumn[i]; |
1253 | | if (!array[lastIndex]) |
1254 | | array[lastIndex] = 1.0e-100; |
1255 | | } |
1256 | | } |
1257 | | } else { |
1258 | | // apply scaling |
1259 | | double scale = model->columnScale()[iBasic]; |
1260 | | for (i = columnStart[iBasic]; |
1261 | | i < columnStart[iBasic] + columnLength[iBasic]; i++) { |
1262 | | int iRow = row[i]; |
1263 | | while (iRow > next) { |
1264 | | lastIndex++; |
1265 | | if (lastIndex == numberOld) |
1266 | | next = matrix_->getNumRows(); |
1267 | | else |
1268 | | next = index[lastIndex]; |
1269 | | } |
1270 | | if (iRow < next) { |
1271 | | array[number] = -elementByColumn[i] * scale * rowScale[iRow]; |
1272 | | index[number++] = iRow; |
1273 | | } else { |
1274 | | assert (iRow == next); |
1275 | | array[lastIndex] -= elementByColumn[i] * scale * rowScale[iRow]; |
1276 | | if (!array[lastIndex]) |
1277 | | array[lastIndex] = 1.0e-100; |
1278 | | } |
1279 | | } |
1280 | | } |
1281 | | rowArray->setNumElements(number); |
1282 | | } |
1283 | | } |
1284 | | } else { |
1285 | | // key slack entering |
1286 | | int iBasic = keyVariable_[gubSlackIn_]; |
1287 | | assert (iBasic < numberColumns); |
1288 | | int number = 0; |
1289 | | const double * rowScale = model->rowScale(); |
1290 | | const int * row = matrix_->getIndices(); |
1291 | | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
1292 | | const int * columnLength = matrix_->getVectorLengths(); |
1293 | | const double * elementByColumn = matrix_->getElements(); |
1294 | | double * array = rowArray->denseVector(); |
1295 | | int * index = rowArray->getIndices(); |
1296 | | CoinBigIndex i; |
1297 | | if (!rowScale) { |
1298 | | for (i = columnStart[iBasic]; |
1299 | | i < columnStart[iBasic] + columnLength[iBasic]; i++) { |
1300 | | int iRow = row[i]; |
1301 | | array[number] = elementByColumn[i]; |
1302 | | index[number++] = iRow; |
1303 | | } |
1304 | | } else { |
1305 | | // apply scaling |
1306 | | double scale = model->columnScale()[iBasic]; |
1307 | | for (i = columnStart[iBasic]; |
1308 | | i < columnStart[iBasic] + columnLength[iBasic]; i++) { |
1309 | | int iRow = row[i]; |
1310 | | array[number] = elementByColumn[i] * scale * rowScale[iRow]; |
1311 | | index[number++] = iRow; |
1312 | | } |
1313 | | } |
1314 | | rowArray->setNumElements(number); |
1315 | | rowArray->setPacked(); |
1316 | | } |
| 1205 | int numberColumns = model->numberColumns(); |
| 1206 | if (iColumn < numberColumns) { |
| 1207 | // Do packed part |
| 1208 | ClpPackedMatrix::unpackPacked(model, rowArray, iColumn); |
| 1209 | int iSet = backward_[iColumn]; |
| 1210 | if (iSet >= 0) { |
| 1211 | // columns are in order |
| 1212 | int iBasic = keyVariable_[iSet]; |
| 1213 | if (iBasic < numberColumns) { |
| 1214 | int number = rowArray->getNumElements(); |
| 1215 | const double *rowScale = model->rowScale(); |
| 1216 | const int *row = matrix_->getIndices(); |
| 1217 | const CoinBigIndex *columnStart = matrix_->getVectorStarts(); |
| 1218 | const int *columnLength = matrix_->getVectorLengths(); |
| 1219 | const double *elementByColumn = matrix_->getElements(); |
| 1220 | double *array = rowArray->denseVector(); |
| 1221 | int *index = rowArray->getIndices(); |
| 1222 | CoinBigIndex i; |
| 1223 | int numberOld = number; |
| 1224 | int lastIndex = 0; |
| 1225 | int next = index[lastIndex]; |
| 1226 | if (!rowScale) { |
| 1227 | for (i = columnStart[iBasic]; |
| 1228 | i < columnStart[iBasic] + columnLength[iBasic]; i++) { |
| 1229 | int iRow = row[i]; |
| 1230 | while (iRow > next) { |
| 1231 | lastIndex++; |
| 1232 | if (lastIndex == numberOld) |
| 1233 | next = matrix_->getNumRows(); |
| 1234 | else |
| 1235 | next = index[lastIndex]; |
| 1236 | } |
| 1237 | if (iRow < next) { |
| 1238 | array[number] = -elementByColumn[i]; |
| 1239 | index[number++] = iRow; |
| 1240 | } else { |
| 1241 | assert(iRow == next); |
| 1242 | array[lastIndex] -= elementByColumn[i]; |
| 1243 | if (!array[lastIndex]) |
| 1244 | array[lastIndex] = 1.0e-100; |
| 1245 | } |
| 1246 | } |
| 1247 | } else { |
| 1248 | // apply scaling |
| 1249 | double scale = model->columnScale()[iBasic]; |
| 1250 | for (i = columnStart[iBasic]; |
| 1251 | i < columnStart[iBasic] + columnLength[iBasic]; i++) { |
| 1252 | int iRow = row[i]; |
| 1253 | while (iRow > next) { |
| 1254 | lastIndex++; |
| 1255 | if (lastIndex == numberOld) |
| 1256 | next = matrix_->getNumRows(); |
| 1257 | else |
| 1258 | next = index[lastIndex]; |
| 1259 | } |
| 1260 | if (iRow < next) { |
| 1261 | array[number] = -elementByColumn[i] * scale * rowScale[iRow]; |
| 1262 | index[number++] = iRow; |
| 1263 | } else { |
| 1264 | assert(iRow == next); |
| 1265 | array[lastIndex] -= elementByColumn[i] * scale * rowScale[iRow]; |
| 1266 | if (!array[lastIndex]) |
| 1267 | array[lastIndex] = 1.0e-100; |
| 1268 | } |
| 1269 | } |
| 1270 | } |
| 1271 | rowArray->setNumElements(number); |
| 1272 | } |
| 1273 | } |
| 1274 | } else { |
| 1275 | // key slack entering |
| 1276 | int iBasic = keyVariable_[gubSlackIn_]; |
| 1277 | assert(iBasic < numberColumns); |
| 1278 | int number = 0; |
| 1279 | const double *rowScale = model->rowScale(); |
| 1280 | const int *row = matrix_->getIndices(); |
| 1281 | const CoinBigIndex *columnStart = matrix_->getVectorStarts(); |
| 1282 | const int *columnLength = matrix_->getVectorLengths(); |
| 1283 | const double *elementByColumn = matrix_->getElements(); |
| 1284 | double *array = rowArray->denseVector(); |
| 1285 | int *index = rowArray->getIndices(); |
| 1286 | CoinBigIndex i; |
| 1287 | if (!rowScale) { |
| 1288 | for (i = columnStart[iBasic]; |
| 1289 | i < columnStart[iBasic] + columnLength[iBasic]; i++) { |
| 1290 | int iRow = row[i]; |
| 1291 | array[number] = elementByColumn[i]; |
| 1292 | index[number++] = iRow; |
| 1293 | } |
| 1294 | } else { |
| 1295 | // apply scaling |
| 1296 | double scale = model->columnScale()[iBasic]; |
| 1297 | for (i = columnStart[iBasic]; |
| 1298 | i < columnStart[iBasic] + columnLength[iBasic]; i++) { |
| 1299 | int iRow = row[i]; |
| 1300 | array[number] = elementByColumn[i] * scale * rowScale[iRow]; |
| 1301 | index[number++] = iRow; |
| 1302 | } |
| 1303 | } |
| 1304 | rowArray->setNumElements(number); |
| 1305 | rowArray->setPacked(); |
| 1306 | } |
1352 | | numberWanted = currentWanted_; |
1353 | | if (numberSets_) { |
1354 | | // Do packed part before gub |
1355 | | int numberColumns = matrix_->getNumCols(); |
1356 | | double ratio = static_cast<double> (firstGub_) / |
1357 | | static_cast<double> (numberColumns); |
1358 | | ClpPackedMatrix::partialPricing(model, startFraction * ratio, |
1359 | | endFraction * ratio, bestSequence, numberWanted); |
1360 | | if (numberWanted || minimumGoodReducedCosts_ < -1) { |
1361 | | // do gub |
1362 | | const double * element = matrix_->getElements(); |
1363 | | const int * row = matrix_->getIndices(); |
1364 | | const CoinBigIndex * startColumn = matrix_->getVectorStarts(); |
1365 | | const int * length = matrix_->getVectorLengths(); |
1366 | | const double * rowScale = model->rowScale(); |
1367 | | const double * columnScale = model->columnScale(); |
1368 | | int iSequence; |
1369 | | CoinBigIndex j; |
1370 | | double tolerance = model->currentDualTolerance(); |
1371 | | double * reducedCost = model->djRegion(); |
1372 | | const double * duals = model->dualRowSolution(); |
1373 | | const double * cost = model->costRegion(); |
1374 | | double bestDj; |
1375 | | int numberColumns = model->numberColumns(); |
1376 | | int numberRows = model->numberRows(); |
1377 | | if (bestSequence >= 0) |
1378 | | bestDj = fabs(this->reducedCost(model, bestSequence)); |
1379 | | else |
1380 | | bestDj = tolerance; |
1381 | | int sequenceOut = model->sequenceOut(); |
1382 | | int saveSequence = bestSequence; |
1383 | | int startG = firstGub_ + static_cast<int> (startFraction * (lastGub_ - firstGub_)); |
1384 | | int endG = firstGub_ + static_cast<int> (endFraction * (lastGub_ - firstGub_)); |
1385 | | endG = CoinMin(lastGub_, endG + 1); |
1386 | | // If nothing found yet can go all the way to end |
1387 | | int endAll = endG; |
1388 | | if (bestSequence < 0 && !startG) |
1389 | | endAll = lastGub_; |
1390 | | int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_; |
1391 | | int minNeg = minimumGoodReducedCosts_ == -1 ? 5 : minimumGoodReducedCosts_; |
1392 | | int nSets = 0; |
1393 | | int iSet = -1; |
1394 | | double djMod = 0.0; |
1395 | | double infeasibilityCost = model->infeasibilityCost(); |
1396 | | if (rowScale) { |
1397 | | double bestDjMod = 0.0; |
1398 | | // scaled |
1399 | | for (iSequence = startG; iSequence < endAll; iSequence++) { |
1400 | | if (numberWanted + minNeg < originalWanted_ && nSets > minSet) { |
1401 | | // give up |
1402 | | numberWanted = 0; |
1403 | | break; |
1404 | | } else if (iSequence == endG && bestSequence >= 0) { |
1405 | | break; |
1406 | | } |
1407 | | if (backward_[iSequence] != iSet) { |
1408 | | // get pi on gub row |
1409 | | iSet = backward_[iSequence]; |
1410 | | if (iSet >= 0) { |
1411 | | nSets++; |
1412 | | int iBasic = keyVariable_[iSet]; |
1413 | | if (iBasic >= numberColumns) { |
1414 | | djMod = - weight(iSet) * infeasibilityCost; |
1415 | | } else { |
1416 | | // get dj without |
1417 | | assert (model->getStatus(iBasic) == ClpSimplex::basic); |
1418 | | djMod = 0.0; |
1419 | | // scaled |
1420 | | for (j = startColumn[iBasic]; |
1421 | | j < startColumn[iBasic] + length[iBasic]; j++) { |
1422 | | int jRow = row[j]; |
1423 | | djMod -= duals[jRow] * element[j] * rowScale[jRow]; |
1424 | | } |
1425 | | // allow for scaling |
1426 | | djMod += cost[iBasic] / columnScale[iBasic]; |
1427 | | // See if gub slack possible - dj is djMod |
1428 | | if (getStatus(iSet) == ClpSimplex::atLowerBound) { |
1429 | | double value = -djMod; |
1430 | | if (value > tolerance) { |
1431 | | numberWanted--; |
1432 | | if (value > bestDj) { |
1433 | | // check flagged variable and correct dj |
1434 | | if (!flagged(iSet)) { |
1435 | | bestDj = value; |
1436 | | bestSequence = numberRows + numberColumns + iSet; |
1437 | | bestDjMod = djMod; |
1438 | | } else { |
1439 | | // just to make sure we don't exit before got something |
1440 | | numberWanted++; |
1441 | | abort(); |
1442 | | } |
1443 | | } |
1444 | | } |
1445 | | } else if (getStatus(iSet) == ClpSimplex::atUpperBound) { |
1446 | | double value = djMod; |
1447 | | if (value > tolerance) { |
1448 | | numberWanted--; |
1449 | | if (value > bestDj) { |
1450 | | // check flagged variable and correct dj |
1451 | | if (!flagged(iSet)) { |
1452 | | bestDj = value; |
1453 | | bestSequence = numberRows + numberColumns + iSet; |
1454 | | bestDjMod = djMod; |
1455 | | } else { |
1456 | | // just to make sure we don't exit before got something |
1457 | | numberWanted++; |
1458 | | abort(); |
1459 | | } |
1460 | | } |
1461 | | } |
1462 | | } |
1463 | | } |
1464 | | } else { |
1465 | | // not in set |
1466 | | djMod = 0.0; |
1467 | | } |
1468 | | } |
1469 | | if (iSequence != sequenceOut) { |
1470 | | double value; |
1471 | | ClpSimplex::Status status = model->getStatus(iSequence); |
| 1339 | numberWanted = currentWanted_; |
| 1340 | if (numberSets_) { |
| 1341 | // Do packed part before gub |
| 1342 | int numberColumns = matrix_->getNumCols(); |
| 1343 | double ratio = static_cast< double >(firstGub_) / static_cast< double >(numberColumns); |
| 1344 | ClpPackedMatrix::partialPricing(model, startFraction * ratio, |
| 1345 | endFraction * ratio, bestSequence, numberWanted); |
| 1346 | if (numberWanted || minimumGoodReducedCosts_ < -1) { |
| 1347 | // do gub |
| 1348 | const double *element = matrix_->getElements(); |
| 1349 | const int *row = matrix_->getIndices(); |
| 1350 | const CoinBigIndex *startColumn = matrix_->getVectorStarts(); |
| 1351 | const int *length = matrix_->getVectorLengths(); |
| 1352 | const double *rowScale = model->rowScale(); |
| 1353 | const double *columnScale = model->columnScale(); |
| 1354 | int iSequence; |
| 1355 | CoinBigIndex j; |
| 1356 | double tolerance = model->currentDualTolerance(); |
| 1357 | double *reducedCost = model->djRegion(); |
| 1358 | const double *duals = model->dualRowSolution(); |
| 1359 | const double *cost = model->costRegion(); |
| 1360 | double bestDj; |
| 1361 | int numberColumns = model->numberColumns(); |
| 1362 | int numberRows = model->numberRows(); |
| 1363 | if (bestSequence >= 0) |
| 1364 | bestDj = fabs(this->reducedCost(model, bestSequence)); |
| 1365 | else |
| 1366 | bestDj = tolerance; |
| 1367 | int sequenceOut = model->sequenceOut(); |
| 1368 | int saveSequence = bestSequence; |
| 1369 | int startG = firstGub_ + static_cast< int >(startFraction * (lastGub_ - firstGub_)); |
| 1370 | int endG = firstGub_ + static_cast< int >(endFraction * (lastGub_ - firstGub_)); |
| 1371 | endG = CoinMin(lastGub_, endG + 1); |
| 1372 | // If nothing found yet can go all the way to end |
| 1373 | int endAll = endG; |
| 1374 | if (bestSequence < 0 && !startG) |
| 1375 | endAll = lastGub_; |
| 1376 | int minSet = minimumObjectsScan_ < 0 ? 5 : minimumObjectsScan_; |
| 1377 | int minNeg = minimumGoodReducedCosts_ == -1 ? 5 : minimumGoodReducedCosts_; |
| 1378 | int nSets = 0; |
| 1379 | int iSet = -1; |
| 1380 | double djMod = 0.0; |
| 1381 | double infeasibilityCost = model->infeasibilityCost(); |
| 1382 | if (rowScale) { |
| 1383 | double bestDjMod = 0.0; |
| 1384 | // scaled |
| 1385 | for (iSequence = startG; iSequence < endAll; iSequence++) { |
| 1386 | if (numberWanted + minNeg < originalWanted_ && nSets > minSet) { |
| 1387 | // give up |
| 1388 | numberWanted = 0; |
| 1389 | break; |
| 1390 | } else if (iSequence == endG && bestSequence >= 0) { |
| 1391 | break; |
| 1392 | } |
| 1393 | if (backward_[iSequence] != iSet) { |
| 1394 | // get pi on gub row |
| 1395 | iSet = backward_[iSequence]; |
| 1396 | if (iSet >= 0) { |
| 1397 | nSets++; |
| 1398 | int iBasic = keyVariable_[iSet]; |
| 1399 | if (iBasic >= numberColumns) { |
| 1400 | djMod = -weight(iSet) * infeasibilityCost; |
| 1401 | } else { |
| 1402 | // get dj without |
| 1403 | assert(model->getStatus(iBasic) == ClpSimplex::basic); |
| 1404 | djMod = 0.0; |
| 1405 | // scaled |
| 1406 | for (j = startColumn[iBasic]; |
| 1407 | j < startColumn[iBasic] + length[iBasic]; j++) { |
| 1408 | int jRow = row[j]; |
| 1409 | djMod -= duals[jRow] * element[j] * rowScale[jRow]; |
| 1410 | } |
| 1411 | // allow for scaling |
| 1412 | djMod += cost[iBasic] / columnScale[iBasic]; |
| 1413 | // See if gub slack possible - dj is djMod |
| 1414 | if (getStatus(iSet) == ClpSimplex::atLowerBound) { |
| 1415 | double value = -djMod; |
| 1416 | if (value > tolerance) { |
| 1417 | numberWanted--; |
| 1418 | if (value > bestDj) { |
| 1419 | // check flagged variable and correct dj |
| 1420 | if (!flagged(iSet)) { |
| 1421 | bestDj = value; |
| 1422 | bestSequence = numberRows + numberColumns + iSet; |
| 1423 | bestDjMod = djMod; |
| 1424 | } else { |
| 1425 | // just to make sure we don't exit before got something |
| 1426 | numberWanted++; |
| 1427 | abort(); |
| 1428 | } |
| 1429 | } |
| 1430 | } |
| 1431 | } else if (getStatus(iSet) == ClpSimplex::atUpperBound) { |
| 1432 | double value = djMod; |
| 1433 | if (value > tolerance) { |
| 1434 | numberWanted--; |
| 1435 | if (value > bestDj) { |
| 1436 | // check flagged variable and correct dj |
| 1437 | if (!flagged(iSet)) { |
| 1438 | bestDj = value; |
| 1439 | bestSequence = numberRows + numberColumns + iSet; |
| 1440 | bestDjMod = djMod; |
| 1441 | } else { |
| 1442 | // just to make sure we don't exit before got something |
| 1443 | numberWanted++; |
| 1444 | abort(); |
| 1445 | } |
| 1446 | } |
| 1447 | } |
| 1448 | } |
| 1449 | } |
| 1450 | } else { |
| 1451 | // not in set |
| 1452 | djMod = 0.0; |
| 1453 | } |
| 1454 | } |
| 1455 | if (iSequence != sequenceOut) { |
| 1456 | double value; |
| 1457 | ClpSimplex::Status status = model->getStatus(iSequence); |
1475 | | case ClpSimplex::basic: |
1476 | | case ClpSimplex::isFixed: |
1477 | | break; |
1478 | | case ClpSimplex::isFree: |
1479 | | case ClpSimplex::superBasic: |
1480 | | value = -djMod; |
1481 | | // scaled |
1482 | | for (j = startColumn[iSequence]; |
1483 | | j < startColumn[iSequence] + length[iSequence]; j++) { |
1484 | | int jRow = row[j]; |
1485 | | value -= duals[jRow] * element[j] * rowScale[jRow]; |
1486 | | } |
1487 | | value = fabs(cost[iSequence] + value * columnScale[iSequence]); |
1488 | | if (value > FREE_ACCEPT * tolerance) { |
1489 | | numberWanted--; |
1490 | | // we are going to bias towards free (but only if reasonable) |
1491 | | value *= FREE_BIAS; |
1492 | | if (value > bestDj) { |
1493 | | // check flagged variable and correct dj |
1494 | | if (!model->flagged(iSequence)) { |
1495 | | bestDj = value; |
1496 | | bestSequence = iSequence; |
1497 | | bestDjMod = djMod; |
1498 | | } else { |
1499 | | // just to make sure we don't exit before got something |
1500 | | numberWanted++; |
1501 | | } |
1502 | | } |
1503 | | } |
1504 | | break; |
1505 | | case ClpSimplex::atUpperBound: |
1506 | | value = -djMod; |
1507 | | // scaled |
1508 | | for (j = startColumn[iSequence]; |
1509 | | j < startColumn[iSequence] + length[iSequence]; j++) { |
1510 | | int jRow = row[j]; |
1511 | | value -= duals[jRow] * element[j] * rowScale[jRow]; |
1512 | | } |
1513 | | value = cost[iSequence] + value * columnScale[iSequence]; |
1514 | | if (value > tolerance) { |
1515 | | numberWanted--; |
1516 | | if (value > bestDj) { |
1517 | | // check flagged variable and correct dj |
1518 | | if (!model->flagged(iSequence)) { |
1519 | | bestDj = value; |
1520 | | bestSequence = iSequence; |
1521 | | bestDjMod = djMod; |
1522 | | } else { |
1523 | | // just to make sure we don't exit before got something |
1524 | | numberWanted++; |
1525 | | } |
1526 | | } |
1527 | | } |
1528 | | break; |
1529 | | case ClpSimplex::atLowerBound: |
1530 | | value = -djMod; |
1531 | | // scaled |
1532 | | for (j = startColumn[iSequence]; |
1533 | | j < startColumn[iSequence] + length[iSequence]; j++) { |
1534 | | int jRow = row[j]; |
1535 | | value -= duals[jRow] * element[j] * rowScale[jRow]; |
1536 | | } |
1537 | | value = -(cost[iSequence] + value * columnScale[iSequence]); |
1538 | | if (value > tolerance) { |
1539 | | numberWanted--; |
1540 | | if (value > bestDj) { |
1541 | | // check flagged variable and correct dj |
1542 | | if (!model->flagged(iSequence)) { |
1543 | | bestDj = value; |
1544 | | bestSequence = iSequence; |
1545 | | bestDjMod = djMod; |
1546 | | } else { |
1547 | | // just to make sure we don't exit before got something |
1548 | | numberWanted++; |
1549 | | } |
1550 | | } |
1551 | | } |
1552 | | break; |
1553 | | } |
1554 | | } |
1555 | | if (!numberWanted) |
1556 | | break; |
| 1461 | case ClpSimplex::basic: |
| 1462 | case ClpSimplex::isFixed: |
| 1463 | break; |
| 1464 | case ClpSimplex::isFree: |
| 1465 | case ClpSimplex::superBasic: |
| 1466 | value = -djMod; |
| 1467 | // scaled |
| 1468 | for (j = startColumn[iSequence]; |
| 1469 | j < startColumn[iSequence] + length[iSequence]; j++) { |
| 1470 | int jRow = row[j]; |
| 1471 | value -= duals[jRow] * element[j] * rowScale[jRow]; |
| 1472 | } |
| 1473 | value = fabs(cost[iSequence] + value * columnScale[iSequence]); |
| 1474 | if (value > FREE_ACCEPT * tolerance) { |
| 1475 | numberWanted--; |
| 1476 | // we are going to bias towards free (but only if reasonable) |
| 1477 | value *= FREE_BIAS; |
| 1478 | if (value > bestDj) { |
| 1479 | // check flagged variable and correct dj |
| 1480 | if (!model->flagged(iSequence)) { |
| 1481 | bestDj = value; |
| 1482 | bestSequence = iSequence; |
| 1483 | bestDjMod = djMod; |
| 1484 | } else { |
| 1485 | // just to make sure we don't exit before got something |
| 1486 | numberWanted++; |
| 1487 | } |
| 1488 | } |
| 1489 | } |
| 1490 | break; |
| 1491 | case ClpSimplex::atUpperBound: |
| 1492 | value = -djMod; |
| 1493 | // scaled |
| 1494 | for (j = startColumn[iSequence]; |
| 1495 | j < startColumn[iSequence] + length[iSequence]; j++) { |
| 1496 | int jRow = row[j]; |
| 1497 | value -= duals[jRow] * element[j] * rowScale[jRow]; |
| 1498 | } |
| 1499 | value = cost[iSequence] + value * columnScale[iSequence]; |
| 1500 | if (value > tolerance) { |
| 1501 | numberWanted--; |
| 1502 | if (value > bestDj) { |
| 1503 | // check flagged variable and correct dj |
| 1504 | if (!model->flagged(iSequence)) { |
| 1505 | bestDj = value; |
| 1506 | bestSequence = iSequence; |
| 1507 | bestDjMod = djMod; |
| 1508 | } else { |
| 1509 | // just to make sure we don't exit before got something |
| 1510 | numberWanted++; |
| 1511 | } |
| 1512 | } |
| 1513 | } |
| 1514 | break; |
| 1515 | case ClpSimplex::atLowerBound: |
| 1516 | value = -djMod; |
| 1517 | // scaled |
| 1518 | for (j = startColumn[iSequence]; |
| 1519 | j < startColumn[iSequence] + length[iSequence]; j++) { |
| 1520 | int jRow = row[j]; |
| 1521 | value -= duals[jRow] * element[j] * rowScale[jRow]; |
| 1522 | } |
| 1523 | value = -(cost[iSequence] + value * columnScale[iSequence]); |
| 1524 | if (value > tolerance) { |
| 1525 | numberWanted--; |
| 1526 | if (value > bestDj) { |
| 1527 | // check flagged variable and correct dj |
| 1528 | if (!model->flagged(iSequence)) { |
| 1529 | bestDj = value; |
| 1530 | bestSequence = iSequence; |
| 1531 | bestDjMod = djMod; |
| 1532 | } else { |
| 1533 | // just to make sure we don't exit before got something |
| 1534 | numberWanted++; |
| 1535 | } |
| 1536 | } |
| 1537 | } |
| 1538 | break; |
| 1539 | } |
| 1540 | } |
| 1541 | if (!numberWanted) |
| 1542 | break; |
| 1543 | } |
| 1544 | if (bestSequence != saveSequence) { |
| 1545 | if (bestSequence < numberRows + numberColumns) { |
| 1546 | // recompute dj |
| 1547 | double value = bestDjMod; |
| 1548 | // scaled |
| 1549 | for (j = startColumn[bestSequence]; |
| 1550 | j < startColumn[bestSequence] + length[bestSequence]; j++) { |
| 1551 | int jRow = row[j]; |
| 1552 | value -= duals[jRow] * element[j] * rowScale[jRow]; |
| 1553 | } |
| 1554 | reducedCost[bestSequence] = cost[bestSequence] + value * columnScale[bestSequence]; |
| 1555 | gubSlackIn_ = -1; |
| 1556 | } else { |
| 1557 | // slack - make last column |
| 1558 | gubSlackIn_ = bestSequence - numberRows - numberColumns; |
| 1559 | bestSequence = numberColumns + 2 * numberRows; |
| 1560 | reducedCost[bestSequence] = bestDjMod; |
| 1561 | model->setStatus(bestSequence, getStatus(gubSlackIn_)); |
| 1562 | if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound) |
| 1563 | model->solutionRegion()[bestSequence] = upper_[gubSlackIn_]; |
| 1564 | else |
| 1565 | model->solutionRegion()[bestSequence] = lower_[gubSlackIn_]; |
| 1566 | model->lowerRegion()[bestSequence] = lower_[gubSlackIn_]; |
| 1567 | model->upperRegion()[bestSequence] = upper_[gubSlackIn_]; |
| 1568 | model->costRegion()[bestSequence] = 0.0; |
| 1569 | } |
| 1570 | savedBestSequence_ = bestSequence; |
| 1571 | savedBestDj_ = reducedCost[savedBestSequence_]; |
| 1572 | } |
| 1573 | } else { |
| 1574 | double bestDjMod = 0.0; |
| 1575 | //printf("iteration %d start %d end %d - wanted %d\n",model->numberIterations(), |
| 1576 | // startG,endG,numberWanted); |
| 1577 | for (iSequence = startG; iSequence < endG; iSequence++) { |
| 1578 | if (numberWanted + minNeg < originalWanted_ && nSets > minSet) { |
| 1579 | // give up |
| 1580 | numberWanted = 0; |
| 1581 | break; |
| 1582 | } else if (iSequence == endG && bestSequence >= 0) { |
| 1583 | break; |
| 1584 | } |
| 1585 | if (backward_[iSequence] != iSet) { |
| 1586 | // get pi on gub row |
| 1587 | iSet = backward_[iSequence]; |
| 1588 | if (iSet >= 0) { |
| 1589 | nSets++; |
| 1590 | int iBasic = keyVariable_[iSet]; |
| 1591 | if (iBasic >= numberColumns) { |
| 1592 | djMod = -weight(iSet) * infeasibilityCost; |
| 1593 | } else { |
| 1594 | // get dj without |
| 1595 | assert(model->getStatus(iBasic) == ClpSimplex::basic); |
| 1596 | djMod = 0.0; |
| 1597 | |
| 1598 | for (j = startColumn[iBasic]; |
| 1599 | j < startColumn[iBasic] + length[iBasic]; j++) { |
| 1600 | int jRow = row[j]; |
| 1601 | djMod -= duals[jRow] * element[j]; |
| 1602 | } |
| 1603 | djMod += cost[iBasic]; |
| 1604 | // See if gub slack possible - dj is djMod |
| 1605 | if (getStatus(iSet) == ClpSimplex::atLowerBound) { |
| 1606 | double value = -djMod; |
| 1607 | if (value > tolerance) { |
| 1608 | numberWanted--; |
| 1609 | if (value > bestDj) { |
| 1610 | // check flagged variable and correct dj |
| 1611 | if (!flagged(iSet)) { |
| 1612 | bestDj = value; |
| 1613 | bestSequence = numberRows + numberColumns + iSet; |
| 1614 | bestDjMod = djMod; |
| 1615 | } else { |
| 1616 | // just to make sure we don't exit before got something |
| 1617 | numberWanted++; |
| 1618 | abort(); |
| 1619 | } |
1664 | | switch(status) { |
1665 | | |
1666 | | case ClpSimplex::basic: |
1667 | | case ClpSimplex::isFixed: |
1668 | | break; |
1669 | | case ClpSimplex::isFree: |
1670 | | case ClpSimplex::superBasic: |
1671 | | value = cost[iSequence] - djMod; |
1672 | | for (j = startColumn[iSequence]; |
1673 | | j < startColumn[iSequence] + length[iSequence]; j++) { |
1674 | | int jRow = row[j]; |
1675 | | value -= duals[jRow] * element[j]; |
1676 | | } |
1677 | | value = fabs(value); |
1678 | | if (value > FREE_ACCEPT * tolerance) { |
1679 | | numberWanted--; |
1680 | | // we are going to bias towards free (but only if reasonable) |
1681 | | value *= FREE_BIAS; |
1682 | | if (value > bestDj) { |
1683 | | // check flagged variable and correct dj |
1684 | | if (!model->flagged(iSequence)) { |
1685 | | bestDj = value; |
1686 | | bestSequence = iSequence; |
1687 | | bestDjMod = djMod; |
1688 | | } else { |
1689 | | // just to make sure we don't exit before got something |
1690 | | numberWanted++; |
1691 | | } |
1692 | | } |
1693 | | } |
1694 | | break; |
1695 | | case ClpSimplex::atUpperBound: |
1696 | | value = cost[iSequence] - djMod; |
1697 | | for (j = startColumn[iSequence]; |
1698 | | j < startColumn[iSequence] + length[iSequence]; j++) { |
1699 | | int jRow = row[j]; |
1700 | | value -= duals[jRow] * element[j]; |
1701 | | } |
1702 | | if (value > tolerance) { |
1703 | | numberWanted--; |
1704 | | if (value > bestDj) { |
1705 | | // check flagged variable and correct dj |
1706 | | if (!model->flagged(iSequence)) { |
1707 | | bestDj = value; |
1708 | | bestSequence = iSequence; |
1709 | | bestDjMod = djMod; |
1710 | | } else { |
1711 | | // just to make sure we don't exit before got something |
1712 | | numberWanted++; |
1713 | | } |
1714 | | } |
1715 | | } |
1716 | | break; |
1717 | | case ClpSimplex::atLowerBound: |
1718 | | value = cost[iSequence] - djMod; |
1719 | | for (j = startColumn[iSequence]; |
1720 | | j < startColumn[iSequence] + length[iSequence]; j++) { |
1721 | | int jRow = row[j]; |
1722 | | value -= duals[jRow] * element[j]; |
1723 | | } |
1724 | | value = -value; |
1725 | | if (value > tolerance) { |
1726 | | numberWanted--; |
1727 | | if (value > bestDj) { |
1728 | | // check flagged variable and correct dj |
1729 | | if (!model->flagged(iSequence)) { |
1730 | | bestDj = value; |
1731 | | bestSequence = iSequence; |
1732 | | bestDjMod = djMod; |
1733 | | } else { |
1734 | | // just to make sure we don't exit before got something |
1735 | | numberWanted++; |
1736 | | } |
1737 | | } |
1738 | | } |
1739 | | break; |
1740 | | } |
1741 | | } |
1742 | | if (!numberWanted) |
1743 | | break; |
1744 | | } |
1745 | | if (bestSequence != saveSequence) { |
1746 | | if (bestSequence < numberRows + numberColumns) { |
1747 | | // recompute dj |
1748 | | double value = cost[bestSequence] - bestDjMod; |
1749 | | for (j = startColumn[bestSequence]; |
1750 | | j < startColumn[bestSequence] + length[bestSequence]; j++) { |
1751 | | int jRow = row[j]; |
1752 | | value -= duals[jRow] * element[j]; |
1753 | | } |
1754 | | //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod); |
1755 | | reducedCost[bestSequence] = value; |
1756 | | gubSlackIn_ = -1; |
1757 | | } else { |
1758 | | // slack - make last column |
1759 | | gubSlackIn_ = bestSequence - numberRows - numberColumns; |
1760 | | bestSequence = numberColumns + 2 * numberRows; |
1761 | | reducedCost[bestSequence] = bestDjMod; |
1762 | | //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod); |
1763 | | model->setStatus(bestSequence, getStatus(gubSlackIn_)); |
1764 | | if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound) |
1765 | | model->solutionRegion()[bestSequence] = upper_[gubSlackIn_]; |
1766 | | else |
1767 | | model->solutionRegion()[bestSequence] = lower_[gubSlackIn_]; |
1768 | | model->lowerRegion()[bestSequence] = lower_[gubSlackIn_]; |
1769 | | model->upperRegion()[bestSequence] = upper_[gubSlackIn_]; |
1770 | | model->costRegion()[bestSequence] = 0.0; |
1771 | | } |
1772 | | } |
1773 | | } |
1774 | | // See if may be finished |
1775 | | if (startG == firstGub_ && bestSequence < 0) |
1776 | | infeasibilityWeight_ = model_->infeasibilityCost(); |
1777 | | else if (bestSequence >= 0) |
1778 | | infeasibilityWeight_ = -1.0; |
1779 | | } |
1780 | | if (numberWanted) { |
1781 | | // Do packed part after gub |
1782 | | double offset = static_cast<double> (lastGub_) / |
1783 | | static_cast<double> (numberColumns); |
1784 | | double ratio = static_cast<double> (numberColumns) / |
1785 | | static_cast<double> (numberColumns) - offset; |
1786 | | double start2 = offset + ratio * startFraction; |
1787 | | double end2 = CoinMin(1.0, offset + ratio * endFraction + 1.0e-6); |
1788 | | ClpPackedMatrix::partialPricing(model, start2, end2, bestSequence, numberWanted); |
1789 | | } |
1790 | | } else { |
1791 | | // no gub |
1792 | | ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted); |
1793 | | } |
1794 | | if (bestSequence >= 0) |
1795 | | infeasibilityWeight_ = -1.0; // not optimal |
1796 | | currentWanted_ = numberWanted; |
| 1652 | case ClpSimplex::basic: |
| 1653 | case ClpSimplex::isFixed: |
| 1654 | break; |
| 1655 | case ClpSimplex::isFree: |
| 1656 | case ClpSimplex::superBasic: |
| 1657 | value = cost[iSequence] - djMod; |
| 1658 | for (j = startColumn[iSequence]; |
| 1659 | j < startColumn[iSequence] + length[iSequence]; j++) { |
| 1660 | int jRow = row[j]; |
| 1661 | value -= duals[jRow] * element[j]; |
| 1662 | } |
| 1663 | value = fabs(value); |
| 1664 | if (value > FREE_ACCEPT * tolerance) { |
| 1665 | numberWanted--; |
| 1666 | // we are going to bias towards free (but only if reasonable) |
| 1667 | value *= FREE_BIAS; |
| 1668 | if (value > bestDj) { |
| 1669 | // check flagged variable and correct dj |
| 1670 | if (!model->flagged(iSequence)) { |
| 1671 | bestDj = value; |
| 1672 | bestSequence = iSequence; |
| 1673 | bestDjMod = djMod; |
| 1674 | } else { |
| 1675 | // just to make sure we don't exit before got something |
| 1676 | numberWanted++; |
| 1677 | } |
| 1678 | } |
| 1679 | } |
| 1680 | break; |
| 1681 | case ClpSimplex::atUpperBound: |
| 1682 | value = cost[iSequence] - djMod; |
| 1683 | for (j = startColumn[iSequence]; |
| 1684 | j < startColumn[iSequence] + length[iSequence]; j++) { |
| 1685 | int jRow = row[j]; |
| 1686 | value -= duals[jRow] * element[j]; |
| 1687 | } |
| 1688 | if (value > tolerance) { |
| 1689 | numberWanted--; |
| 1690 | if (value > bestDj) { |
| 1691 | // check flagged variable and correct dj |
| 1692 | if (!model->flagged(iSequence)) { |
| 1693 | bestDj = value; |
| 1694 | bestSequence = iSequence; |
| 1695 | bestDjMod = djMod; |
| 1696 | } else { |
| 1697 | // just to make sure we don't exit before got something |
| 1698 | numberWanted++; |
| 1699 | } |
| 1700 | } |
| 1701 | } |
| 1702 | break; |
| 1703 | case ClpSimplex::atLowerBound: |
| 1704 | value = cost[iSequence] - djMod; |
| 1705 | for (j = startColumn[iSequence]; |
| 1706 | j < startColumn[iSequence] + length[iSequence]; j++) { |
| 1707 | int jRow = row[j]; |
| 1708 | value -= duals[jRow] * element[j]; |
| 1709 | } |
| 1710 | value = -value; |
| 1711 | if (value > tolerance) { |
| 1712 | numberWanted--; |
| 1713 | if (value > bestDj) { |
| 1714 | // check flagged variable and correct dj |
| 1715 | if (!model->flagged(iSequence)) { |
| 1716 | bestDj = value; |
| 1717 | bestSequence = iSequence; |
| 1718 | bestDjMod = djMod; |
| 1719 | } else { |
| 1720 | // just to make sure we don't exit before got something |
| 1721 | numberWanted++; |
| 1722 | } |
| 1723 | } |
| 1724 | } |
| 1725 | break; |
| 1726 | } |
| 1727 | } |
| 1728 | if (!numberWanted) |
| 1729 | break; |
| 1730 | } |
| 1731 | if (bestSequence != saveSequence) { |
| 1732 | if (bestSequence < numberRows + numberColumns) { |
| 1733 | // recompute dj |
| 1734 | double value = cost[bestSequence] - bestDjMod; |
| 1735 | for (j = startColumn[bestSequence]; |
| 1736 | j < startColumn[bestSequence] + length[bestSequence]; j++) { |
| 1737 | int jRow = row[j]; |
| 1738 | value -= duals[jRow] * element[j]; |
| 1739 | } |
| 1740 | //printf("price struct %d - dj %g gubpi %g\n",bestSequence,value,bestDjMod); |
| 1741 | reducedCost[bestSequence] = value; |
| 1742 | gubSlackIn_ = -1; |
| 1743 | } else { |
| 1744 | // slack - make last column |
| 1745 | gubSlackIn_ = bestSequence - numberRows - numberColumns; |
| 1746 | bestSequence = numberColumns + 2 * numberRows; |
| 1747 | reducedCost[bestSequence] = bestDjMod; |
| 1748 | //printf("price slack %d - gubpi %g\n",gubSlackIn_,bestDjMod); |
| 1749 | model->setStatus(bestSequence, getStatus(gubSlackIn_)); |
| 1750 | if (getStatus(gubSlackIn_) == ClpSimplex::atUpperBound) |
| 1751 | model->solutionRegion()[bestSequence] = upper_[gubSlackIn_]; |
| 1752 | else |
| 1753 | model->solutionRegion()[bestSequence] = lower_[gubSlackIn_]; |
| 1754 | model->lowerRegion()[bestSequence] = lower_[gubSlackIn_]; |
| 1755 | model->upperRegion()[bestSequence] = upper_[gubSlackIn_]; |
| 1756 | model->costRegion()[bestSequence] = 0.0; |
| 1757 | } |
| 1758 | } |
| 1759 | } |
| 1760 | // See if may be finished |
| 1761 | if (startG == firstGub_ && bestSequence < 0) |
| 1762 | infeasibilityWeight_ = model_->infeasibilityCost(); |
| 1763 | else if (bestSequence >= 0) |
| 1764 | infeasibilityWeight_ = -1.0; |
| 1765 | } |
| 1766 | if (numberWanted) { |
| 1767 | // Do packed part after gub |
| 1768 | double offset = static_cast< double >(lastGub_) / static_cast< double >(numberColumns); |
| 1769 | double ratio = static_cast< double >(numberColumns) / static_cast< double >(numberColumns) - offset; |
| 1770 | double start2 = offset + ratio * startFraction; |
| 1771 | double end2 = CoinMin(1.0, offset + ratio * endFraction + 1.0e-6); |
| 1772 | ClpPackedMatrix::partialPricing(model, start2, end2, bestSequence, numberWanted); |
| 1773 | } |
| 1774 | } else { |
| 1775 | // no gub |
| 1776 | ClpPackedMatrix::partialPricing(model, startFraction, endFraction, bestSequence, numberWanted); |
| 1777 | } |
| 1778 | if (bestSequence >= 0) |
| 1779 | infeasibilityWeight_ = -1.0; // not optimal |
| 1780 | currentWanted_ = numberWanted; |
1804 | | // I think we only need to bother about sets with two in basis or incoming set |
1805 | | int number = update->getNumElements(); |
1806 | | double * array = update->denseVector(); |
1807 | | int * index = update->getIndices(); |
1808 | | int i; |
1809 | | assert (!number || update->packedMode()); |
1810 | | int * pivotVariable = model->pivotVariable(); |
1811 | | int numberRows = model->numberRows(); |
1812 | | int numberColumns = model->numberColumns(); |
1813 | | int numberTotal = numberRows + numberColumns; |
1814 | | int sequenceIn = model->sequenceIn(); |
1815 | | int returnCode = 0; |
1816 | | int iSetIn; |
1817 | | if (sequenceIn < numberColumns) { |
1818 | | iSetIn = backward_[sequenceIn]; |
1819 | | gubSlackIn_ = -1; // in case set |
1820 | | } else if (sequenceIn < numberRows + numberColumns) { |
1821 | | iSetIn = -1; |
1822 | | gubSlackIn_ = -1; // in case set |
1823 | | } else { |
1824 | | iSetIn = gubSlackIn_; |
1825 | | } |
1826 | | double * lower = model->lowerRegion(); |
1827 | | double * upper = model->upperRegion(); |
1828 | | double * cost = model->costRegion(); |
1829 | | double * solution = model->solutionRegion(); |
1830 | | int number2 = number; |
1831 | | if (!mode) { |
1832 | | double primalTolerance = model->primalTolerance(); |
1833 | | double infeasibilityCost = model->infeasibilityCost(); |
1834 | | // extend |
1835 | | saveNumber_ = number; |
1836 | | for (i = 0; i < number; i++) { |
1837 | | int iRow = index[i]; |
1838 | | int iPivot = pivotVariable[iRow]; |
1839 | | if (iPivot < numberColumns) { |
1840 | | int iSet = backward_[iPivot]; |
1841 | | if (iSet >= 0) { |
1842 | | // two (or more) in set |
1843 | | int iIndex = toIndex_[iSet]; |
1844 | | double otherValue = array[i]; |
1845 | | double value; |
1846 | | if (iIndex < 0) { |
1847 | | toIndex_[iSet] = number2; |
1848 | | int iNew = number2 - number; |
1849 | | fromIndex_[number2-number] = iSet; |
1850 | | iIndex = number2; |
1851 | | index[number2] = numberRows + iNew; |
1852 | | // do key stuff |
1853 | | int iKey = keyVariable_[iSet]; |
1854 | | if (iKey < numberColumns) { |
1855 | | // Save current cost of key |
1856 | | changeCost_[number2-number] = cost[iKey]; |
1857 | | if (iSet != iSetIn) |
1858 | | value = 0.0; |
1859 | | else if (iSetIn != gubSlackIn_) |
1860 | | value = 1.0; |
1861 | | else |
1862 | | value = -1.0; |
1863 | | pivotVariable[numberRows+iNew] = iKey; |
1864 | | // Do I need to recompute? |
1865 | | double sol; |
1866 | | assert (getStatus(iSet) != ClpSimplex::basic); |
1867 | | if (getStatus(iSet) == ClpSimplex::atLowerBound) |
1868 | | sol = lower_[iSet]; |
1869 | | else |
1870 | | sol = upper_[iSet]; |
1871 | | if ((gubType_ & 8) != 0) { |
1872 | | int iColumn = next_[iKey]; |
1873 | | // sum all non-key variables |
1874 | | while(iColumn >= 0) { |
1875 | | sol -= solution[iColumn]; |
1876 | | iColumn = next_[iColumn]; |
1877 | | } |
1878 | | } else { |
1879 | | int stop = -(iKey + 1); |
1880 | | int iColumn = next_[iKey]; |
1881 | | // sum all non-key variables |
1882 | | while(iColumn != stop) { |
1883 | | if (iColumn < 0) |
1884 | | iColumn = -iColumn - 1; |
1885 | | sol -= solution[iColumn]; |
1886 | | iColumn = next_[iColumn]; |
1887 | | } |
1888 | | } |
1889 | | solution[iKey] = sol; |
1890 | | if (model->algorithm() > 0) |
1891 | | model->nonLinearCost()->setOne(iKey, sol); |
1892 | | //assert (fabs(sol-solution[iKey])<1.0e-3); |
1893 | | } else { |
1894 | | // gub slack is basic |
1895 | | // Save current cost of key |
1896 | | changeCost_[number2-number] = -weight(iSet) * infeasibilityCost; |
1897 | | otherValue = - otherValue; //allow for - sign on slack |
1898 | | if (iSet != iSetIn) |
1899 | | value = 0.0; |
1900 | | else |
1901 | | value = -1.0; |
1902 | | pivotVariable[numberRows+iNew] = iNew + numberTotal; |
1903 | | model->djRegion()[iNew+numberTotal] = 0.0; |
1904 | | double sol = 0.0; |
1905 | | if ((gubType_ & 8) != 0) { |
1906 | | int iColumn = next_[iKey]; |
1907 | | // sum all non-key variables |
1908 | | while(iColumn >= 0) { |
1909 | | sol += solution[iColumn]; |
1910 | | iColumn = next_[iColumn]; |
1911 | | } |
1912 | | } else { |
1913 | | int stop = -(iKey + 1); |
1914 | | int iColumn = next_[iKey]; |
1915 | | // sum all non-key variables |
1916 | | while(iColumn != stop) { |
1917 | | if (iColumn < 0) |
1918 | | iColumn = -iColumn - 1; |
1919 | | sol += solution[iColumn]; |
1920 | | iColumn = next_[iColumn]; |
1921 | | } |
1922 | | } |
1923 | | solution[iNew+numberTotal] = sol; |
1924 | | // and do cost in nonLinearCost |
1925 | | if (model->algorithm() > 0) |
1926 | | model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSet], upper_[iSet]); |
1927 | | if (sol > upper_[iSet] + primalTolerance) { |
1928 | | setAbove(iSet); |
1929 | | lower[iNew+numberTotal] = upper_[iSet]; |
1930 | | upper[iNew+numberTotal] = COIN_DBL_MAX; |
1931 | | } else if (sol < lower_[iSet] - primalTolerance) { |
1932 | | setBelow(iSet); |
1933 | | lower[iNew+numberTotal] = -COIN_DBL_MAX; |
1934 | | upper[iNew+numberTotal] = lower_[iSet]; |
1935 | | } else { |
1936 | | setFeasible(iSet); |
1937 | | lower[iNew+numberTotal] = lower_[iSet]; |
1938 | | upper[iNew+numberTotal] = upper_[iSet]; |
1939 | | } |
1940 | | cost[iNew+numberTotal] = weight(iSet) * infeasibilityCost; |
1941 | | } |
1942 | | number2++; |
1943 | | } else { |
1944 | | value = array[iIndex]; |
1945 | | int iKey = keyVariable_[iSet]; |
1946 | | if (iKey >= numberColumns) |
1947 | | otherValue = - otherValue; //allow for - sign on slack |
1948 | | } |
1949 | | value -= otherValue; |
1950 | | array[iIndex] = value; |
1951 | | } |
1952 | | } |
1953 | | } |
1954 | | if (iSetIn >= 0 && toIndex_[iSetIn] < 0) { |
1955 | | // Do incoming |
1956 | | update->setPacked(); // just in case no elements |
1957 | | toIndex_[iSetIn] = number2; |
1958 | | int iNew = number2 - number; |
1959 | | fromIndex_[number2-number] = iSetIn; |
1960 | | // Save current cost of key |
1961 | | double currentCost; |
1962 | | int key = keyVariable_[iSetIn]; |
1963 | | if (key < numberColumns) |
1964 | | currentCost = cost[key]; |
1965 | | else |
1966 | | currentCost = -weight(iSetIn) * infeasibilityCost; |
1967 | | changeCost_[number2-number] = currentCost; |
1968 | | index[number2] = numberRows + iNew; |
1969 | | // do key stuff |
1970 | | int iKey = keyVariable_[iSetIn]; |
1971 | | if (iKey < numberColumns) { |
1972 | | if (gubSlackIn_ < 0) |
1973 | | array[number2] = 1.0; |
1974 | | else |
1975 | | array[number2] = -1.0; |
1976 | | pivotVariable[numberRows+iNew] = iKey; |
1977 | | // Do I need to recompute? |
1978 | | double sol; |
1979 | | assert (getStatus(iSetIn) != ClpSimplex::basic); |
1980 | | if (getStatus(iSetIn) == ClpSimplex::atLowerBound) |
1981 | | sol = lower_[iSetIn]; |
1982 | | else |
1983 | | sol = upper_[iSetIn]; |
1984 | | if ((gubType_ & 8) != 0) { |
1985 | | int iColumn = next_[iKey]; |
1986 | | // sum all non-key variables |
1987 | | while(iColumn >= 0) { |
1988 | | sol -= solution[iColumn]; |
1989 | | iColumn = next_[iColumn]; |
1990 | | } |
1991 | | } else { |
1992 | | // bounds exist - sum over all except key |
1993 | | int stop = -(iKey + 1); |
1994 | | int iColumn = next_[iKey]; |
1995 | | // sum all non-key variables |
1996 | | while(iColumn != stop) { |
1997 | | if (iColumn < 0) |
1998 | | iColumn = -iColumn - 1; |
1999 | | sol -= solution[iColumn]; |
2000 | | iColumn = next_[iColumn]; |
2001 | | } |
2002 | | } |
2003 | | solution[iKey] = sol; |
2004 | | if (model->algorithm() > 0) |
2005 | | model->nonLinearCost()->setOne(iKey, sol); |
2006 | | //assert (fabs(sol-solution[iKey])<1.0e-3); |
2007 | | } else { |
2008 | | // gub slack is basic |
2009 | | array[number2] = -1.0; |
2010 | | pivotVariable[numberRows+iNew] = iNew + numberTotal; |
2011 | | model->djRegion()[iNew+numberTotal] = 0.0; |
2012 | | double sol = 0.0; |
2013 | | if ((gubType_ & 8) != 0) { |
2014 | | int iColumn = next_[iKey]; |
2015 | | // sum all non-key variables |
2016 | | while(iColumn >= 0) { |
2017 | | sol += solution[iColumn]; |
2018 | | iColumn = next_[iColumn]; |
2019 | | } |
2020 | | } else { |
2021 | | // bounds exist - sum over all except key |
2022 | | int stop = -(iKey + 1); |
2023 | | int iColumn = next_[iKey]; |
2024 | | // sum all non-key variables |
2025 | | while(iColumn != stop) { |
2026 | | if (iColumn < 0) |
2027 | | iColumn = -iColumn - 1; |
2028 | | sol += solution[iColumn]; |
2029 | | iColumn = next_[iColumn]; |
2030 | | } |
2031 | | } |
2032 | | solution[iNew+numberTotal] = sol; |
2033 | | // and do cost in nonLinearCost |
2034 | | if (model->algorithm() > 0) |
2035 | | model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSetIn], upper_[iSetIn]); |
2036 | | if (sol > upper_[iSetIn] + primalTolerance) { |
2037 | | setAbove(iSetIn); |
2038 | | lower[iNew+numberTotal] = upper_[iSetIn]; |
2039 | | upper[iNew+numberTotal] = COIN_DBL_MAX; |
2040 | | } else if (sol < lower_[iSetIn] - primalTolerance) { |
2041 | | setBelow(iSetIn); |
2042 | | lower[iNew+numberTotal] = -COIN_DBL_MAX; |
2043 | | upper[iNew+numberTotal] = lower_[iSetIn]; |
2044 | | } else { |
2045 | | setFeasible(iSetIn); |
2046 | | lower[iNew+numberTotal] = lower_[iSetIn]; |
2047 | | upper[iNew+numberTotal] = upper_[iSetIn]; |
2048 | | } |
2049 | | cost[iNew+numberTotal] = weight(iSetIn) * infeasibilityCost; |
2050 | | } |
2051 | | number2++; |
2052 | | } |
2053 | | // mark end |
2054 | | fromIndex_[number2-number] = -1; |
2055 | | returnCode = number2 - number; |
2056 | | // make sure lower_ upper_ adjusted |
2057 | | synchronize(model, 9); |
2058 | | } else { |
2059 | | // take off? |
2060 | | if (number > saveNumber_) { |
2061 | | // clear |
2062 | | double theta = model->theta(); |
2063 | | double * solution = model->solutionRegion(); |
2064 | | for (i = saveNumber_; i < number; i++) { |
2065 | | int iRow = index[i]; |
2066 | | int iColumn = pivotVariable[iRow]; |
| 1787 | // I think we only need to bother about sets with two in basis or incoming set |
| 1788 | int number = update->getNumElements(); |
| 1789 | double *array = update->denseVector(); |
| 1790 | int *index = update->getIndices(); |
| 1791 | int i; |
| 1792 | assert(!number || update->packedMode()); |
| 1793 | int *pivotVariable = model->pivotVariable(); |
| 1794 | int numberRows = model->numberRows(); |
| 1795 | int numberColumns = model->numberColumns(); |
| 1796 | int numberTotal = numberRows + numberColumns; |
| 1797 | int sequenceIn = model->sequenceIn(); |
| 1798 | int returnCode = 0; |
| 1799 | int iSetIn; |
| 1800 | if (sequenceIn < numberColumns) { |
| 1801 | iSetIn = backward_[sequenceIn]; |
| 1802 | gubSlackIn_ = -1; // in case set |
| 1803 | } else if (sequenceIn < numberRows + numberColumns) { |
| 1804 | iSetIn = -1; |
| 1805 | gubSlackIn_ = -1; // in case set |
| 1806 | } else { |
| 1807 | iSetIn = gubSlackIn_; |
| 1808 | } |
| 1809 | double *lower = model->lowerRegion(); |
| 1810 | double *upper = model->upperRegion(); |
| 1811 | double *cost = model->costRegion(); |
| 1812 | double *solution = model->solutionRegion(); |
| 1813 | int number2 = number; |
| 1814 | if (!mode) { |
| 1815 | double primalTolerance = model->primalTolerance(); |
| 1816 | double infeasibilityCost = model->infeasibilityCost(); |
| 1817 | // extend |
| 1818 | saveNumber_ = number; |
| 1819 | for (i = 0; i < number; i++) { |
| 1820 | int iRow = index[i]; |
| 1821 | int iPivot = pivotVariable[iRow]; |
| 1822 | if (iPivot < numberColumns) { |
| 1823 | int iSet = backward_[iPivot]; |
| 1824 | if (iSet >= 0) { |
| 1825 | // two (or more) in set |
| 1826 | int iIndex = toIndex_[iSet]; |
| 1827 | double otherValue = array[i]; |
| 1828 | double value; |
| 1829 | if (iIndex < 0) { |
| 1830 | toIndex_[iSet] = number2; |
| 1831 | int iNew = number2 - number; |
| 1832 | fromIndex_[number2 - number] = iSet; |
| 1833 | iIndex = number2; |
| 1834 | index[number2] = numberRows + iNew; |
| 1835 | // do key stuff |
| 1836 | int iKey = keyVariable_[iSet]; |
| 1837 | if (iKey < numberColumns) { |
| 1838 | // Save current cost of key |
| 1839 | changeCost_[number2 - number] = cost[iKey]; |
| 1840 | if (iSet != iSetIn) |
| 1841 | value = 0.0; |
| 1842 | else if (iSetIn != gubSlackIn_) |
| 1843 | value = 1.0; |
| 1844 | else |
| 1845 | value = -1.0; |
| 1846 | pivotVariable[numberRows + iNew] = iKey; |
| 1847 | // Do I need to recompute? |
| 1848 | double sol; |
| 1849 | assert(getStatus(iSet) != ClpSimplex::basic); |
| 1850 | if (getStatus(iSet) == ClpSimplex::atLowerBound) |
| 1851 | sol = lower_[iSet]; |
| 1852 | else |
| 1853 | sol = upper_[iSet]; |
| 1854 | if ((gubType_ & 8) != 0) { |
| 1855 | int iColumn = next_[iKey]; |
| 1856 | // sum all non-key variables |
| 1857 | while (iColumn >= 0) { |
| 1858 | sol -= solution[iColumn]; |
| 1859 | iColumn = next_[iColumn]; |
| 1860 | } |
| 1861 | } else { |
| 1862 | int stop = -(iKey + 1); |
| 1863 | int iColumn = next_[iKey]; |
| 1864 | // sum all non-key variables |
| 1865 | while (iColumn != stop) { |
| 1866 | if (iColumn < 0) |
| 1867 | iColumn = -iColumn - 1; |
| 1868 | sol -= solution[iColumn]; |
| 1869 | iColumn = next_[iColumn]; |
| 1870 | } |
| 1871 | } |
| 1872 | solution[iKey] = sol; |
| 1873 | if (model->algorithm() > 0) |
| 1874 | model->nonLinearCost()->setOne(iKey, sol); |
| 1875 | //assert (fabs(sol-solution[iKey])<1.0e-3); |
| 1876 | } else { |
| 1877 | // gub slack is basic |
| 1878 | // Save current cost of key |
| 1879 | changeCost_[number2 - number] = -weight(iSet) * infeasibilityCost; |
| 1880 | otherValue = -otherValue; //allow for - sign on slack |
| 1881 | if (iSet != iSetIn) |
| 1882 | value = 0.0; |
| 1883 | else |
| 1884 | value = -1.0; |
| 1885 | pivotVariable[numberRows + iNew] = iNew + numberTotal; |
| 1886 | model->djRegion()[iNew + numberTotal] = 0.0; |
| 1887 | double sol = 0.0; |
| 1888 | if ((gubType_ & 8) != 0) { |
| 1889 | int iColumn = next_[iKey]; |
| 1890 | // sum all non-key variables |
| 1891 | while (iColumn >= 0) { |
| 1892 | sol += solution[iColumn]; |
| 1893 | iColumn = next_[iColumn]; |
| 1894 | } |
| 1895 | } else { |
| 1896 | int stop = -(iKey + 1); |
| 1897 | int iColumn = next_[iKey]; |
| 1898 | // sum all non-key variables |
| 1899 | while (iColumn != stop) { |
| 1900 | if (iColumn < 0) |
| 1901 | iColumn = -iColumn - 1; |
| 1902 | sol += solution[iColumn]; |
| 1903 | iColumn = next_[iColumn]; |
| 1904 | } |
| 1905 | } |
| 1906 | solution[iNew + numberTotal] = sol; |
| 1907 | // and do cost in nonLinearCost |
| 1908 | if (model->algorithm() > 0) |
| 1909 | model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSet], upper_[iSet]); |
| 1910 | if (sol > upper_[iSet] + primalTolerance) { |
| 1911 | setAbove(iSet); |
| 1912 | lower[iNew + numberTotal] = upper_[iSet]; |
| 1913 | upper[iNew + numberTotal] = COIN_DBL_MAX; |
| 1914 | } else if (sol < lower_[iSet] - primalTolerance) { |
| 1915 | setBelow(iSet); |
| 1916 | lower[iNew + numberTotal] = -COIN_DBL_MAX; |
| 1917 | upper[iNew + numberTotal] = lower_[iSet]; |
| 1918 | } else { |
| 1919 | setFeasible(iSet); |
| 1920 | lower[iNew + numberTotal] = lower_[iSet]; |
| 1921 | upper[iNew + numberTotal] = upper_[iSet]; |
| 1922 | } |
| 1923 | cost[iNew + numberTotal] = weight(iSet) * infeasibilityCost; |
| 1924 | } |
| 1925 | number2++; |
| 1926 | } else { |
| 1927 | value = array[iIndex]; |
| 1928 | int iKey = keyVariable_[iSet]; |
| 1929 | if (iKey >= numberColumns) |
| 1930 | otherValue = -otherValue; //allow for - sign on slack |
| 1931 | } |
| 1932 | value -= otherValue; |
| 1933 | array[iIndex] = value; |
| 1934 | } |
| 1935 | } |
| 1936 | } |
| 1937 | if (iSetIn >= 0 && toIndex_[iSetIn] < 0) { |
| 1938 | // Do incoming |
| 1939 | update->setPacked(); // just in case no elements |
| 1940 | toIndex_[iSetIn] = number2; |
| 1941 | int iNew = number2 - number; |
| 1942 | fromIndex_[number2 - number] = iSetIn; |
| 1943 | // Save current cost of key |
| 1944 | double currentCost; |
| 1945 | int key = keyVariable_[iSetIn]; |
| 1946 | if (key < numberColumns) |
| 1947 | currentCost = cost[key]; |
| 1948 | else |
| 1949 | currentCost = -weight(iSetIn) * infeasibilityCost; |
| 1950 | changeCost_[number2 - number] = currentCost; |
| 1951 | index[number2] = numberRows + iNew; |
| 1952 | // do key stuff |
| 1953 | int iKey = keyVariable_[iSetIn]; |
| 1954 | if (iKey < numberColumns) { |
| 1955 | if (gubSlackIn_ < 0) |
| 1956 | array[number2] = 1.0; |
| 1957 | else |
| 1958 | array[number2] = -1.0; |
| 1959 | pivotVariable[numberRows + iNew] = iKey; |
| 1960 | // Do I need to recompute? |
| 1961 | double sol; |
| 1962 | assert(getStatus(iSetIn) != ClpSimplex::basic); |
| 1963 | if (getStatus(iSetIn) == ClpSimplex::atLowerBound) |
| 1964 | sol = lower_[iSetIn]; |
| 1965 | else |
| 1966 | sol = upper_[iSetIn]; |
| 1967 | if ((gubType_ & 8) != 0) { |
| 1968 | int iColumn = next_[iKey]; |
| 1969 | // sum all non-key variables |
| 1970 | while (iColumn >= 0) { |
| 1971 | sol -= solution[iColumn]; |
| 1972 | iColumn = next_[iColumn]; |
| 1973 | } |
| 1974 | } else { |
| 1975 | // bounds exist - sum over all except key |
| 1976 | int stop = -(iKey + 1); |
| 1977 | int iColumn = next_[iKey]; |
| 1978 | // sum all non-key variables |
| 1979 | while (iColumn != stop) { |
| 1980 | if (iColumn < 0) |
| 1981 | iColumn = -iColumn - 1; |
| 1982 | sol -= solution[iColumn]; |
| 1983 | iColumn = next_[iColumn]; |
| 1984 | } |
| 1985 | } |
| 1986 | solution[iKey] = sol; |
| 1987 | if (model->algorithm() > 0) |
| 1988 | model->nonLinearCost()->setOne(iKey, sol); |
| 1989 | //assert (fabs(sol-solution[iKey])<1.0e-3); |
| 1990 | } else { |
| 1991 | // gub slack is basic |
| 1992 | array[number2] = -1.0; |
| 1993 | pivotVariable[numberRows + iNew] = iNew + numberTotal; |
| 1994 | model->djRegion()[iNew + numberTotal] = 0.0; |
| 1995 | double sol = 0.0; |
| 1996 | if ((gubType_ & 8) != 0) { |
| 1997 | int iColumn = next_[iKey]; |
| 1998 | // sum all non-key variables |
| 1999 | while (iColumn >= 0) { |
| 2000 | sol += solution[iColumn]; |
| 2001 | iColumn = next_[iColumn]; |
| 2002 | } |
| 2003 | } else { |
| 2004 | // bounds exist - sum over all except key |
| 2005 | int stop = -(iKey + 1); |
| 2006 | int iColumn = next_[iKey]; |
| 2007 | // sum all non-key variables |
| 2008 | while (iColumn != stop) { |
| 2009 | if (iColumn < 0) |
| 2010 | iColumn = -iColumn - 1; |
| 2011 | sol += solution[iColumn]; |
| 2012 | iColumn = next_[iColumn]; |
| 2013 | } |
| 2014 | } |
| 2015 | solution[iNew + numberTotal] = sol; |
| 2016 | // and do cost in nonLinearCost |
| 2017 | if (model->algorithm() > 0) |
| 2018 | model->nonLinearCost()->setOne(iNew + numberTotal, sol, lower_[iSetIn], upper_[iSetIn]); |
| 2019 | if (sol > upper_[iSetIn] + primalTolerance) { |
| 2020 | setAbove(iSetIn); |
| 2021 | lower[iNew + numberTotal] = upper_[iSetIn]; |
| 2022 | upper[iNew + numberTotal] = COIN_DBL_MAX; |
| 2023 | } else if (sol < lower_[iSetIn] - primalTolerance) { |
| 2024 | setBelow(iSetIn); |
| 2025 | lower[iNew + numberTotal] = -COIN_DBL_MAX; |
| 2026 | upper[iNew + numberTotal] = lower_[iSetIn]; |
| 2027 | } else { |
| 2028 | setFeasible(iSetIn); |
| 2029 | lower[iNew + numberTotal] = lower_[iSetIn]; |
| 2030 | upper[iNew + numberTotal] = upper_[iSetIn]; |
| 2031 | } |
| 2032 | cost[iNew + numberTotal] = weight(iSetIn) * infeasibilityCost; |
| 2033 | } |
| 2034 | number2++; |
| 2035 | } |
| 2036 | // mark end |
| 2037 | fromIndex_[number2 - number] = -1; |
| 2038 | returnCode = number2 - number; |
| 2039 | // make sure lower_ upper_ adjusted |
| 2040 | synchronize(model, 9); |
| 2041 | } else { |
| 2042 | // take off? |
| 2043 | if (number > saveNumber_) { |
| 2044 | // clear |
| 2045 | double theta = model->theta(); |
| 2046 | double *solution = model->solutionRegion(); |
| 2047 | for (i = saveNumber_; i < number; i++) { |
| 2048 | int iRow = index[i]; |
| 2049 | int iColumn = pivotVariable[iRow]; |
2099 | | int numberColumns = model->numberColumns(); |
2100 | | switch (mode) { |
2101 | | // If key variable then slot in gub rhs so will get correct contribution |
2102 | | case 0: { |
2103 | | int i; |
2104 | | double * solution = model->solutionRegion(); |
2105 | | ClpSimplex::Status iStatus; |
2106 | | for (i = 0; i < numberSets_; i++) { |
2107 | | int iColumn = keyVariable_[i]; |
2108 | | if (iColumn < numberColumns) { |
2109 | | // key is structural - where is slack |
2110 | | iStatus = getStatus(i); |
2111 | | assert (iStatus != ClpSimplex::basic); |
2112 | | if (iStatus == ClpSimplex::atLowerBound) |
2113 | | solution[iColumn] = lower_[i]; |
2114 | | else |
2115 | | solution[iColumn] = upper_[i]; |
2116 | | } |
2117 | | } |
2118 | | } |
2119 | | break; |
2120 | | // Compute values of key variables |
2121 | | case 1: { |
2122 | | int i; |
2123 | | double * solution = model->solutionRegion(); |
2124 | | //const int * columnLength = matrix_->getVectorLengths(); |
2125 | | //const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
2126 | | //const int * row = matrix_->getIndices(); |
2127 | | //const double * elementByColumn = matrix_->getElements(); |
2128 | | //int * pivotVariable = model->pivotVariable(); |
2129 | | sumPrimalInfeasibilities_ = 0.0; |
2130 | | numberPrimalInfeasibilities_ = 0; |
2131 | | double primalTolerance = model->primalTolerance(); |
2132 | | double relaxedTolerance = primalTolerance; |
2133 | | // we can't really trust infeasibilities if there is primal error |
2134 | | double error = CoinMin(1.0e-2, model->largestPrimalError()); |
2135 | | // allow tolerance at least slightly bigger than standard |
2136 | | relaxedTolerance = relaxedTolerance + error; |
2137 | | // but we will be using difference |
2138 | | relaxedTolerance -= primalTolerance; |
2139 | | sumOfRelaxedPrimalInfeasibilities_ = 0.0; |
2140 | | for (i = 0; i < numberSets_; i++) { // Could just be over basics (esp if no bounds) |
2141 | | int kColumn = keyVariable_[i]; |
2142 | | double value = 0.0; |
2143 | | if ((gubType_ & 8) != 0) { |
2144 | | int iColumn = next_[kColumn]; |
2145 | | // sum all non-key variables |
2146 | | while(iColumn >= 0) { |
2147 | | value += solution[iColumn]; |
2148 | | iColumn = next_[iColumn]; |
2149 | | } |
2150 | | } else { |
2151 | | // bounds exist - sum over all except key |
2152 | | int stop = -(kColumn + 1); |
2153 | | int iColumn = next_[kColumn]; |
2154 | | // sum all non-key variables |
2155 | | while(iColumn != stop) { |
2156 | | if (iColumn < 0) |
2157 | | iColumn = -iColumn - 1; |
2158 | | value += solution[iColumn]; |
2159 | | iColumn = next_[iColumn]; |
2160 | | } |
2161 | | } |
2162 | | if (kColumn < numberColumns) { |
2163 | | // make sure key is basic - so will be skipped in values pass |
2164 | | model->setStatus(kColumn, ClpSimplex::basic); |
2165 | | // feasibility will be done later |
2166 | | assert (getStatus(i) != ClpSimplex::basic); |
2167 | | if (getStatus(i) == ClpSimplex::atUpperBound) |
2168 | | solution[kColumn] = upper_[i] - value; |
2169 | | else |
2170 | | solution[kColumn] = lower_[i] - value; |
2171 | | //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]); |
2172 | | } else { |
2173 | | // slack is key |
2174 | | assert (getStatus(i) == ClpSimplex::basic); |
2175 | | double infeasibility = 0.0; |
2176 | | if (value > upper_[i] + primalTolerance) { |
2177 | | infeasibility = value - upper_[i] - primalTolerance; |
2178 | | setAbove(i); |
2179 | | } else if (value < lower_[i] - primalTolerance) { |
2180 | | infeasibility = lower_[i] - value - primalTolerance ; |
2181 | | setBelow(i); |
2182 | | } else { |
2183 | | setFeasible(i); |
2184 | | } |
2185 | | //printf("Value of key slack for set %d is %g\n",i,value); |
2186 | | if (infeasibility > 0.0) { |
2187 | | sumPrimalInfeasibilities_ += infeasibility; |
2188 | | if (infeasibility > relaxedTolerance) |
2189 | | sumOfRelaxedPrimalInfeasibilities_ += infeasibility; |
2190 | | numberPrimalInfeasibilities_ ++; |
2191 | | } |
2192 | | } |
2193 | | } |
2194 | | } |
2195 | | break; |
2196 | | // Report on infeasibilities of key variables |
2197 | | case 2: { |
2198 | | model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities() + |
2199 | | sumPrimalInfeasibilities_); |
2200 | | model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities() + |
2201 | | numberPrimalInfeasibilities_); |
2202 | | model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities() + |
2203 | | sumOfRelaxedPrimalInfeasibilities_); |
2204 | | } |
2205 | | break; |
2206 | | } |
| 2081 | int numberColumns = model->numberColumns(); |
| 2082 | switch (mode) { |
| 2083 | // If key variable then slot in gub rhs so will get correct contribution |
| 2084 | case 0: { |
| 2085 | int i; |
| 2086 | double *solution = model->solutionRegion(); |
| 2087 | ClpSimplex::Status iStatus; |
| 2088 | for (i = 0; i < numberSets_; i++) { |
| 2089 | int iColumn = keyVariable_[i]; |
| 2090 | if (iColumn < numberColumns) { |
| 2091 | // key is structural - where is slack |
| 2092 | iStatus = getStatus(i); |
| 2093 | assert(iStatus != ClpSimplex::basic); |
| 2094 | if (iStatus == ClpSimplex::atLowerBound) |
| 2095 | solution[iColumn] = lower_[i]; |
| 2096 | else |
| 2097 | solution[iColumn] = upper_[i]; |
| 2098 | } |
| 2099 | } |
| 2100 | } break; |
| 2101 | // Compute values of key variables |
| 2102 | case 1: { |
| 2103 | int i; |
| 2104 | double *solution = model->solutionRegion(); |
| 2105 | //const int * columnLength = matrix_->getVectorLengths(); |
| 2106 | //const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
| 2107 | //const int * row = matrix_->getIndices(); |
| 2108 | //const double * elementByColumn = matrix_->getElements(); |
| 2109 | //int * pivotVariable = model->pivotVariable(); |
| 2110 | sumPrimalInfeasibilities_ = 0.0; |
| 2111 | numberPrimalInfeasibilities_ = 0; |
| 2112 | double primalTolerance = model->primalTolerance(); |
| 2113 | double relaxedTolerance = primalTolerance; |
| 2114 | // we can't really trust infeasibilities if there is primal error |
| 2115 | double error = CoinMin(1.0e-2, model->largestPrimalError()); |
| 2116 | // allow tolerance at least slightly bigger than standard |
| 2117 | relaxedTolerance = relaxedTolerance + error; |
| 2118 | // but we will be using difference |
| 2119 | relaxedTolerance -= primalTolerance; |
| 2120 | sumOfRelaxedPrimalInfeasibilities_ = 0.0; |
| 2121 | for (i = 0; i < numberSets_; i++) { // Could just be over basics (esp if no bounds) |
| 2122 | int kColumn = keyVariable_[i]; |
| 2123 | double value = 0.0; |
| 2124 | if ((gubType_ & 8) != 0) { |
| 2125 | int iColumn = next_[kColumn]; |
| 2126 | // sum all non-key variables |
| 2127 | while (iColumn >= 0) { |
| 2128 | value += solution[iColumn]; |
| 2129 | iColumn = next_[iColumn]; |
| 2130 | } |
| 2131 | } else { |
| 2132 | // bounds exist - sum over all except key |
| 2133 | int stop = -(kColumn + 1); |
| 2134 | int iColumn = next_[kColumn]; |
| 2135 | // sum all non-key variables |
| 2136 | while (iColumn != stop) { |
| 2137 | if (iColumn < 0) |
| 2138 | iColumn = -iColumn - 1; |
| 2139 | value += solution[iColumn]; |
| 2140 | iColumn = next_[iColumn]; |
| 2141 | } |
| 2142 | } |
| 2143 | if (kColumn < numberColumns) { |
| 2144 | // make sure key is basic - so will be skipped in values pass |
| 2145 | model->setStatus(kColumn, ClpSimplex::basic); |
| 2146 | // feasibility will be done later |
| 2147 | assert(getStatus(i) != ClpSimplex::basic); |
| 2148 | if (getStatus(i) == ClpSimplex::atUpperBound) |
| 2149 | solution[kColumn] = upper_[i] - value; |
| 2150 | else |
| 2151 | solution[kColumn] = lower_[i] - value; |
| 2152 | //printf("Value of key structural %d for set %d is %g\n",kColumn,i,solution[kColumn]); |
| 2153 | } else { |
| 2154 | // slack is key |
| 2155 | assert(getStatus(i) == ClpSimplex::basic); |
| 2156 | double infeasibility = 0.0; |
| 2157 | if (value > upper_[i] + primalTolerance) { |
| 2158 | infeasibility = value - upper_[i] - primalTolerance; |
| 2159 | setAbove(i); |
| 2160 | } else if (value < lower_[i] - primalTolerance) { |
| 2161 | infeasibility = lower_[i] - value - primalTolerance; |
| 2162 | setBelow(i); |
| 2163 | } else { |
| 2164 | setFeasible(i); |
| 2165 | } |
| 2166 | //printf("Value of key slack for set %d is %g\n",i,value); |
| 2167 | if (infeasibility > 0.0) { |
| 2168 | sumPrimalInfeasibilities_ += infeasibility; |
| 2169 | if (infeasibility > relaxedTolerance) |
| 2170 | sumOfRelaxedPrimalInfeasibilities_ += infeasibility; |
| 2171 | numberPrimalInfeasibilities_++; |
| 2172 | } |
| 2173 | } |
| 2174 | } |
| 2175 | } break; |
| 2176 | // Report on infeasibilities of key variables |
| 2177 | case 2: { |
| 2178 | model->setSumPrimalInfeasibilities(model->sumPrimalInfeasibilities() + sumPrimalInfeasibilities_); |
| 2179 | model->setNumberPrimalInfeasibilities(model->numberPrimalInfeasibilities() + numberPrimalInfeasibilities_); |
| 2180 | model->setSumOfRelaxedPrimalInfeasibilities(model->sumOfRelaxedPrimalInfeasibilities() + sumOfRelaxedPrimalInfeasibilities_); |
| 2181 | } break; |
| 2182 | } |
2218 | | switch (mode) { |
2219 | | // modify costs before transposeUpdate |
2220 | | case 0: { |
2221 | | int i; |
2222 | | double * cost = model->costRegion(); |
2223 | | // not dual values yet |
2224 | | //assert (!other); |
2225 | | //double * work = array->denseVector(); |
2226 | | double infeasibilityCost = model->infeasibilityCost(); |
2227 | | int * pivotVariable = model->pivotVariable(); |
2228 | | int numberRows = model->numberRows(); |
2229 | | int numberColumns = model->numberColumns(); |
2230 | | for (i = 0; i < numberRows; i++) { |
2231 | | int iPivot = pivotVariable[i]; |
2232 | | if (iPivot < numberColumns) { |
2233 | | int iSet = backward_[iPivot]; |
2234 | | if (iSet >= 0) { |
2235 | | int kColumn = keyVariable_[iSet]; |
2236 | | double costValue; |
2237 | | if (kColumn < numberColumns) { |
2238 | | // structural has cost |
2239 | | costValue = cost[kColumn]; |
2240 | | } else { |
2241 | | // slack is key |
2242 | | assert (getStatus(iSet) == ClpSimplex::basic); |
2243 | | // negative as -1.0 for slack |
2244 | | costValue = -weight(iSet) * infeasibilityCost; |
2245 | | } |
2246 | | array->add(i, -costValue); // was work[i]-costValue |
2247 | | } |
2248 | | } |
2249 | | } |
2250 | | } |
2251 | | break; |
2252 | | // create duals for key variables (without check on dual infeasible) |
2253 | | case 1: { |
2254 | | // If key slack then dual 0.0 (if feasible) |
2255 | | // dj for key is zero so that defines dual on set |
2256 | | int i; |
2257 | | double * dj = model->djRegion(); |
2258 | | int numberColumns = model->numberColumns(); |
2259 | | double infeasibilityCost = model->infeasibilityCost(); |
2260 | | for (i = 0; i < numberSets_; i++) { |
2261 | | int kColumn = keyVariable_[i]; |
2262 | | if (kColumn < numberColumns) { |
2263 | | // dj without set |
2264 | | double value = dj[kColumn]; |
2265 | | // Now subtract out from all |
2266 | | dj[kColumn] = 0.0; |
2267 | | int iColumn = next_[kColumn]; |
2268 | | // modify all non-key variables |
2269 | | while(iColumn >= 0) { |
2270 | | dj[iColumn] -= value; |
2271 | | iColumn = next_[iColumn]; |
2272 | | } |
2273 | | } else { |
2274 | | // slack key - may not be feasible |
2275 | | assert (getStatus(i) == ClpSimplex::basic); |
2276 | | // negative as -1.0 for slack |
2277 | | double value = -weight(i) * infeasibilityCost; |
2278 | | if (value) { |
2279 | | int iColumn = next_[kColumn]; |
2280 | | // modify all non-key variables basic |
2281 | | while(iColumn >= 0) { |
2282 | | dj[iColumn] -= value; |
2283 | | iColumn = next_[iColumn]; |
2284 | | } |
2285 | | } |
2286 | | } |
2287 | | } |
2288 | | } |
2289 | | break; |
2290 | | // as 1 but check slacks and compute djs |
2291 | | case 2: { |
2292 | | // If key slack then dual 0.0 |
2293 | | // If not then slack could be dual infeasible |
2294 | | // dj for key is zero so that defines dual on set |
2295 | | int i; |
2296 | | // make sure fromIndex will not confuse pricing |
2297 | | fromIndex_[0] = -1; |
| 2193 | switch (mode) { |
| 2194 | // modify costs before transposeUpdate |
| 2195 | case 0: { |
| 2196 | int i; |
| 2197 | double *cost = model->costRegion(); |
| 2198 | // not dual values yet |
| 2199 | //assert (!other); |
| 2200 | //double * work = array->denseVector(); |
| 2201 | double infeasibilityCost = model->infeasibilityCost(); |
| 2202 | int *pivotVariable = model->pivotVariable(); |
| 2203 | int numberRows = model->numberRows(); |
| 2204 | int numberColumns = model->numberColumns(); |
| 2205 | for (i = 0; i < numberRows; i++) { |
| 2206 | int iPivot = pivotVariable[i]; |
| 2207 | if (iPivot < numberColumns) { |
| 2208 | int iSet = backward_[iPivot]; |
| 2209 | if (iSet >= 0) { |
| 2210 | int kColumn = keyVariable_[iSet]; |
| 2211 | double costValue; |
| 2212 | if (kColumn < numberColumns) { |
| 2213 | // structural has cost |
| 2214 | costValue = cost[kColumn]; |
| 2215 | } else { |
| 2216 | // slack is key |
| 2217 | assert(getStatus(iSet) == ClpSimplex::basic); |
| 2218 | // negative as -1.0 for slack |
| 2219 | costValue = -weight(iSet) * infeasibilityCost; |
| 2220 | } |
| 2221 | array->add(i, -costValue); // was work[i]-costValue |
| 2222 | } |
| 2223 | } |
| 2224 | } |
| 2225 | } break; |
| 2226 | // create duals for key variables (without check on dual infeasible) |
| 2227 | case 1: { |
| 2228 | // If key slack then dual 0.0 (if feasible) |
| 2229 | // dj for key is zero so that defines dual on set |
| 2230 | int i; |
| 2231 | double *dj = model->djRegion(); |
| 2232 | int numberColumns = model->numberColumns(); |
| 2233 | double infeasibilityCost = model->infeasibilityCost(); |
| 2234 | for (i = 0; i < numberSets_; i++) { |
| 2235 | int kColumn = keyVariable_[i]; |
| 2236 | if (kColumn < numberColumns) { |
| 2237 | // dj without set |
| 2238 | double value = dj[kColumn]; |
| 2239 | // Now subtract out from all |
| 2240 | dj[kColumn] = 0.0; |
| 2241 | int iColumn = next_[kColumn]; |
| 2242 | // modify all non-key variables |
| 2243 | while (iColumn >= 0) { |
| 2244 | dj[iColumn] -= value; |
| 2245 | iColumn = next_[iColumn]; |
| 2246 | } |
| 2247 | } else { |
| 2248 | // slack key - may not be feasible |
| 2249 | assert(getStatus(i) == ClpSimplex::basic); |
| 2250 | // negative as -1.0 for slack |
| 2251 | double value = -weight(i) * infeasibilityCost; |
| 2252 | if (value) { |
| 2253 | int iColumn = next_[kColumn]; |
| 2254 | // modify all non-key variables basic |
| 2255 | while (iColumn >= 0) { |
| 2256 | dj[iColumn] -= value; |
| 2257 | iColumn = next_[iColumn]; |
| 2258 | } |
| 2259 | } |
| 2260 | } |
| 2261 | } |
| 2262 | } break; |
| 2263 | // as 1 but check slacks and compute djs |
| 2264 | case 2: { |
| 2265 | // If key slack then dual 0.0 |
| 2266 | // If not then slack could be dual infeasible |
| 2267 | // dj for key is zero so that defines dual on set |
| 2268 | int i; |
| 2269 | // make sure fromIndex will not confuse pricing |
| 2270 | fromIndex_[0] = -1; |
| 2271 | possiblePivotKey_ = -1; |
| 2272 | // Create array |
| 2273 | int numberColumns = model->numberColumns(); |
| 2274 | int *pivotVariable = model->pivotVariable(); |
| 2275 | int numberRows = model->numberRows(); |
| 2276 | for (i = 0; i < numberRows; i++) { |
| 2277 | int iPivot = pivotVariable[i]; |
| 2278 | if (iPivot < numberColumns) |
| 2279 | backToPivotRow_[iPivot] = i; |
| 2280 | } |
| 2281 | if (noCheck_ >= 0) { |
| 2282 | if (infeasibilityWeight_ != model->infeasibilityCost()) { |
| 2283 | // don't bother checking |
| 2284 | sumDualInfeasibilities_ = 100.0; |
| 2285 | numberDualInfeasibilities_ = 1; |
| 2286 | sumOfRelaxedDualInfeasibilities_ = 100.0; |
| 2287 | return; |
| 2288 | } |
| 2289 | } |
| 2290 | double *dj = model->djRegion(); |
| 2291 | double *dual = model->dualRowSolution(); |
| 2292 | double *cost = model->costRegion(); |
| 2293 | ClpSimplex::Status iStatus; |
| 2294 | const int *columnLength = matrix_->getVectorLengths(); |
| 2295 | const CoinBigIndex *columnStart = matrix_->getVectorStarts(); |
| 2296 | const int *row = matrix_->getIndices(); |
| 2297 | const double *elementByColumn = matrix_->getElements(); |
| 2298 | double infeasibilityCost = model->infeasibilityCost(); |
| 2299 | sumDualInfeasibilities_ = 0.0; |
| 2300 | numberDualInfeasibilities_ = 0; |
| 2301 | double dualTolerance = model->dualTolerance(); |
| 2302 | double relaxedTolerance = dualTolerance; |
| 2303 | // we can't really trust infeasibilities if there is dual error |
| 2304 | double error = CoinMin(1.0e-2, model->largestDualError()); |
| 2305 | // allow tolerance at least slightly bigger than standard |
| 2306 | relaxedTolerance = relaxedTolerance + error; |
| 2307 | // but we will be using difference |
| 2308 | relaxedTolerance -= dualTolerance; |
| 2309 | sumOfRelaxedDualInfeasibilities_ = 0.0; |
| 2310 | for (i = 0; i < numberSets_; i++) { |
| 2311 | int kColumn = keyVariable_[i]; |
| 2312 | if (kColumn < numberColumns) { |
| 2313 | // dj without set |
| 2314 | double value = cost[kColumn]; |
| 2315 | for (CoinBigIndex j = columnStart[kColumn]; |
| 2316 | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
| 2317 | int iRow = row[j]; |
| 2318 | value -= dual[iRow] * elementByColumn[j]; |
| 2319 | } |
| 2320 | // Now subtract out from all |
| 2321 | dj[kColumn] -= value; |
| 2322 | int stop = -(kColumn + 1); |
| 2323 | kColumn = next_[kColumn]; |
| 2324 | while (kColumn != stop) { |
| 2325 | if (kColumn < 0) |
| 2326 | kColumn = -kColumn - 1; |
| 2327 | double djValue = dj[kColumn] - value; |
| 2328 | dj[kColumn] = djValue; |
| 2329 | ; |
| 2330 | double infeasibility = 0.0; |
| 2331 | iStatus = model->getStatus(kColumn); |
| 2332 | if (iStatus == ClpSimplex::atLowerBound) { |
| 2333 | if (djValue < -dualTolerance) |
| 2334 | infeasibility = -djValue - dualTolerance; |
| 2335 | } else if (iStatus == ClpSimplex::atUpperBound) { |
| 2336 | // at upper bound |
| 2337 | if (djValue > dualTolerance) |
| 2338 | infeasibility = djValue - dualTolerance; |
| 2339 | } |
| 2340 | if (infeasibility > 0.0) { |
| 2341 | sumDualInfeasibilities_ += infeasibility; |
| 2342 | if (infeasibility > relaxedTolerance) |
| 2343 | sumOfRelaxedDualInfeasibilities_ += infeasibility; |
| 2344 | numberDualInfeasibilities_++; |
| 2345 | } |
| 2346 | kColumn = next_[kColumn]; |
| 2347 | } |
| 2348 | // check slack |
| 2349 | iStatus = getStatus(i); |
| 2350 | assert(iStatus != ClpSimplex::basic); |
| 2351 | double infeasibility = 0.0; |
| 2352 | // dj of slack is -(-1.0)value |
| 2353 | if (iStatus == ClpSimplex::atLowerBound) { |
| 2354 | if (value < -dualTolerance) |
| 2355 | infeasibility = -value - dualTolerance; |
| 2356 | } else if (iStatus == ClpSimplex::atUpperBound) { |
| 2357 | // at upper bound |
| 2358 | if (value > dualTolerance) |
| 2359 | infeasibility = value - dualTolerance; |
| 2360 | } |
| 2361 | if (infeasibility > 0.0) { |
| 2362 | sumDualInfeasibilities_ += infeasibility; |
| 2363 | if (infeasibility > relaxedTolerance) |
| 2364 | sumOfRelaxedDualInfeasibilities_ += infeasibility; |
| 2365 | numberDualInfeasibilities_++; |
| 2366 | } |
| 2367 | } else { |
| 2368 | // slack key - may not be feasible |
| 2369 | assert(getStatus(i) == ClpSimplex::basic); |
| 2370 | // negative as -1.0 for slack |
| 2371 | double value = -weight(i) * infeasibilityCost; |
| 2372 | if (value) { |
| 2373 | // Now subtract out from all |
| 2374 | int kColumn = i + numberColumns; |
| 2375 | int stop = -(kColumn + 1); |
| 2376 | kColumn = next_[kColumn]; |
| 2377 | while (kColumn != stop) { |
| 2378 | if (kColumn < 0) |
| 2379 | kColumn = -kColumn - 1; |
| 2380 | double djValue = dj[kColumn] - value; |
| 2381 | dj[kColumn] = djValue; |
| 2382 | ; |
| 2383 | double infeasibility = 0.0; |
| 2384 | iStatus = model->getStatus(kColumn); |
| 2385 | if (iStatus == ClpSimplex::atLowerBound) { |
| 2386 | if (djValue < -dualTolerance) |
| 2387 | infeasibility = -djValue - dualTolerance; |
| 2388 | } else if (iStatus == ClpSimplex::atUpperBound) { |
| 2389 | // at upper bound |
| 2390 | if (djValue > dualTolerance) |
| 2391 | infeasibility = djValue - dualTolerance; |
| 2392 | } |
| 2393 | if (infeasibility > 0.0) { |
| 2394 | sumDualInfeasibilities_ += infeasibility; |
| 2395 | if (infeasibility > relaxedTolerance) |
| 2396 | sumOfRelaxedDualInfeasibilities_ += infeasibility; |
| 2397 | numberDualInfeasibilities_++; |
| 2398 | } |
| 2399 | kColumn = next_[kColumn]; |
| 2400 | } |
| 2401 | } |
| 2402 | } |
| 2403 | } |
| 2404 | // and get statistics for column generation |
| 2405 | synchronize(model, 4); |
| 2406 | infeasibilityWeight_ = -1.0; |
| 2407 | } break; |
| 2408 | // Report on infeasibilities of key variables |
| 2409 | case 3: { |
| 2410 | model->setSumDualInfeasibilities(model->sumDualInfeasibilities() + sumDualInfeasibilities_); |
| 2411 | model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() + numberDualInfeasibilities_); |
| 2412 | model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() + sumOfRelaxedDualInfeasibilities_); |
| 2413 | } break; |
| 2414 | // modify costs before transposeUpdate for partial pricing |
| 2415 | case 4: { |
| 2416 | // First compute new costs etc for interesting gubs |
| 2417 | int iLook = 0; |
| 2418 | int iSet = fromIndex_[0]; |
| 2419 | double primalTolerance = model->primalTolerance(); |
| 2420 | const double *cost = model->costRegion(); |
| 2421 | double *solution = model->solutionRegion(); |
| 2422 | double infeasibilityCost = model->infeasibilityCost(); |
| 2423 | int numberColumns = model->numberColumns(); |
| 2424 | int numberChanged = 0; |
| 2425 | int *pivotVariable = model->pivotVariable(); |
| 2426 | while (iSet >= 0) { |
| 2427 | int key = keyVariable_[iSet]; |
| 2428 | double value = 0.0; |
| 2429 | // sum over all except key |
| 2430 | if ((gubType_ & 8) != 0) { |
| 2431 | int iColumn = next_[key]; |
| 2432 | // sum all non-key variables |
| 2433 | while (iColumn >= 0) { |
| 2434 | value += solution[iColumn]; |
| 2435 | iColumn = next_[iColumn]; |
| 2436 | } |
| 2437 | } else { |
| 2438 | // bounds exist - sum over all except key |
| 2439 | int stop = -(key + 1); |
| 2440 | int iColumn = next_[key]; |
| 2441 | // sum all non-key variables |
| 2442 | while (iColumn != stop) { |
| 2443 | if (iColumn < 0) |
| 2444 | iColumn = -iColumn - 1; |
| 2445 | value += solution[iColumn]; |
| 2446 | iColumn = next_[iColumn]; |
| 2447 | } |
| 2448 | } |
| 2449 | double costChange; |
| 2450 | double oldCost = changeCost_[iLook]; |
| 2451 | if (key < numberColumns) { |
| 2452 | assert(getStatus(iSet) != ClpSimplex::basic); |
| 2453 | double sol; |
| 2454 | if (getStatus(iSet) == ClpSimplex::atUpperBound) |
| 2455 | sol = upper_[iSet] - value; |
| 2456 | else |
| 2457 | sol = lower_[iSet] - value; |
| 2458 | solution[key] = sol; |
| 2459 | // fix up cost |
| 2460 | model->nonLinearCost()->setOne(key, sol); |
| 2461 | #ifdef CLP_DEBUG_PRINT |
| 2462 | printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n", key, iSet, sol, |
| 2463 | cost[key], oldCost); |
| 2464 | #endif |
| 2465 | costChange = cost[key] - oldCost; |
| 2466 | } else { |
| 2467 | // slack is key |
| 2468 | if (value > upper_[iSet] + primalTolerance) { |
| 2469 | setAbove(iSet); |
| 2470 | } else if (value < lower_[iSet] - primalTolerance) { |
| 2471 | setBelow(iSet); |
| 2472 | } else { |
| 2473 | setFeasible(iSet); |
| 2474 | } |
| 2475 | // negative as -1.0 for slack |
| 2476 | costChange = -weight(iSet) * infeasibilityCost - oldCost; |
| 2477 | #ifdef CLP_DEBUG_PRINT |
| 2478 | printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n", iSet, value, |
| 2479 | weight(iSet) * infeasibilityCost, oldCost); |
| 2480 | #endif |
| 2481 | } |
| 2482 | if (costChange) { |
| 2483 | fromIndex_[numberChanged] = iSet; |
| 2484 | toIndex_[iSet] = numberChanged; |
| 2485 | changeCost_[numberChanged++] = costChange; |
| 2486 | } |
| 2487 | iSet = fromIndex_[++iLook]; |
| 2488 | } |
| 2489 | if (numberChanged || possiblePivotKey_ >= 0) { |
| 2490 | // first do those in list already |
| 2491 | int number = array->getNumElements(); |
| 2492 | array->setPacked(); |
| 2493 | int i; |
| 2494 | double *work = array->denseVector(); |
| 2495 | int *which = array->getIndices(); |
| 2496 | for (i = 0; i < number; i++) { |
| 2497 | int iRow = which[i]; |
| 2498 | int iPivot = pivotVariable[iRow]; |
| 2499 | if (iPivot < numberColumns) { |
| 2500 | int iSet = backward_[iPivot]; |
| 2501 | if (iSet >= 0 && toIndex_[iSet] >= 0) { |
| 2502 | double newValue = work[i] + changeCost_[toIndex_[iSet]]; |
| 2503 | if (!newValue) |
| 2504 | newValue = 1.0e-100; |
| 2505 | work[i] = newValue; |
| 2506 | // mark as done |
| 2507 | backward_[iPivot] = -1; |
| 2508 | } |
| 2509 | } |
| 2510 | if (possiblePivotKey_ == iRow) { |
| 2511 | double newValue = work[i] - model->dualIn(); |
| 2512 | if (!newValue) |
| 2513 | newValue = 1.0e-100; |
| 2514 | work[i] = newValue; |
2299 | | // Create array |
2300 | | int numberColumns = model->numberColumns(); |
2301 | | int * pivotVariable = model->pivotVariable(); |
2302 | | int numberRows = model->numberRows(); |
2303 | | for (i = 0; i < numberRows; i++) { |
2304 | | int iPivot = pivotVariable[i]; |
2305 | | if (iPivot < numberColumns) |
2306 | | backToPivotRow_[iPivot] = i; |
2307 | | } |
2308 | | if (noCheck_ >= 0) { |
2309 | | if (infeasibilityWeight_ != model->infeasibilityCost()) { |
2310 | | // don't bother checking |
2311 | | sumDualInfeasibilities_ = 100.0; |
2312 | | numberDualInfeasibilities_ = 1; |
2313 | | sumOfRelaxedDualInfeasibilities_ = 100.0; |
2314 | | return; |
2315 | | } |
2316 | | } |
2317 | | double * dj = model->djRegion(); |
2318 | | double * dual = model->dualRowSolution(); |
2319 | | double * cost = model->costRegion(); |
2320 | | ClpSimplex::Status iStatus; |
2321 | | const int * columnLength = matrix_->getVectorLengths(); |
2322 | | const CoinBigIndex * columnStart = matrix_->getVectorStarts(); |
2323 | | const int * row = matrix_->getIndices(); |
2324 | | const double * elementByColumn = matrix_->getElements(); |
2325 | | double infeasibilityCost = model->infeasibilityCost(); |
2326 | | sumDualInfeasibilities_ = 0.0; |
2327 | | numberDualInfeasibilities_ = 0; |
2328 | | double dualTolerance = model->dualTolerance(); |
2329 | | double relaxedTolerance = dualTolerance; |
2330 | | // we can't really trust infeasibilities if there is dual error |
2331 | | double error = CoinMin(1.0e-2, model->largestDualError()); |
2332 | | // allow tolerance at least slightly bigger than standard |
2333 | | relaxedTolerance = relaxedTolerance + error; |
2334 | | // but we will be using difference |
2335 | | relaxedTolerance -= dualTolerance; |
2336 | | sumOfRelaxedDualInfeasibilities_ = 0.0; |
2337 | | for (i = 0; i < numberSets_; i++) { |
2338 | | int kColumn = keyVariable_[i]; |
2339 | | if (kColumn < numberColumns) { |
2340 | | // dj without set |
2341 | | double value = cost[kColumn]; |
2342 | | for (CoinBigIndex j = columnStart[kColumn]; |
2343 | | j < columnStart[kColumn] + columnLength[kColumn]; j++) { |
2344 | | int iRow = row[j]; |
2345 | | value -= dual[iRow] * elementByColumn[j]; |
2346 | | } |
2347 | | // Now subtract out from all |
2348 | | dj[kColumn] -= value; |
2349 | | int stop = -(kColumn + 1); |
2350 | | kColumn = next_[kColumn]; |
2351 | | while (kColumn != stop) { |
2352 | | if (kColumn < 0) |
2353 | | kColumn = -kColumn - 1; |
2354 | | double djValue = dj[kColumn] - value; |
2355 | | dj[kColumn] = djValue;; |
2356 | | double infeasibility = 0.0; |
2357 | | iStatus = model->getStatus(kColumn); |
2358 | | if (iStatus == ClpSimplex::atLowerBound) { |
2359 | | if (djValue < -dualTolerance) |
2360 | | infeasibility = -djValue - dualTolerance; |
2361 | | } else if (iStatus == ClpSimplex::atUpperBound) { |
2362 | | // at upper bound |
2363 | | if (djValue > dualTolerance) |
2364 | | infeasibility = djValue - dualTolerance; |
2365 | | } |
2366 | | if (infeasibility > 0.0) { |
2367 | | sumDualInfeasibilities_ += infeasibility; |
2368 | | if (infeasibility > relaxedTolerance) |
2369 | | sumOfRelaxedDualInfeasibilities_ += infeasibility; |
2370 | | numberDualInfeasibilities_ ++; |
2371 | | } |
2372 | | kColumn = next_[kColumn]; |
2373 | | } |
2374 | | // check slack |
2375 | | iStatus = getStatus(i); |
2376 | | assert (iStatus != ClpSimplex::basic); |
2377 | | double infeasibility = 0.0; |
2378 | | // dj of slack is -(-1.0)value |
2379 | | if (iStatus == ClpSimplex::atLowerBound) { |
2380 | | if (value < -dualTolerance) |
2381 | | infeasibility = -value - dualTolerance; |
2382 | | } else if (iStatus == ClpSimplex::atUpperBound) { |
2383 | | // at upper bound |
2384 | | if (value > dualTolerance) |
2385 | | infeasibility = value - dualTolerance; |
2386 | | } |
2387 | | if (infeasibility > 0.0) { |
2388 | | sumDualInfeasibilities_ += infeasibility; |
2389 | | if (infeasibility > relaxedTolerance) |
2390 | | sumOfRelaxedDualInfeasibilities_ += infeasibility; |
2391 | | numberDualInfeasibilities_ ++; |
2392 | | } |
2393 | | } else { |
2394 | | // slack key - may not be feasible |
2395 | | assert (getStatus(i) == ClpSimplex::basic); |
2396 | | // negative as -1.0 for slack |
2397 | | double value = -weight(i) * infeasibilityCost; |
2398 | | if (value) { |
2399 | | // Now subtract out from all |
2400 | | int kColumn = i + numberColumns; |
2401 | | int stop = -(kColumn + 1); |
2402 | | kColumn = next_[kColumn]; |
2403 | | while (kColumn != stop) { |
2404 | | if (kColumn < 0) |
2405 | | kColumn = -kColumn - 1; |
2406 | | double djValue = dj[kColumn] - value; |
2407 | | dj[kColumn] = djValue;; |
2408 | | double infeasibility = 0.0; |
2409 | | iStatus = model->getStatus(kColumn); |
2410 | | if (iStatus == ClpSimplex::atLowerBound) { |
2411 | | if (djValue < -dualTolerance) |
2412 | | infeasibility = -djValue - dualTolerance; |
2413 | | } else if (iStatus == ClpSimplex::atUpperBound) { |
2414 | | // at upper bound |
2415 | | if (djValue > dualTolerance) |
2416 | | infeasibility = djValue - dualTolerance; |
2417 | | } |
2418 | | if (infeasibility > 0.0) { |
2419 | | sumDualInfeasibilities_ += infeasibility; |
2420 | | if (infeasibility > relaxedTolerance) |
2421 | | sumOfRelaxedDualInfeasibilities_ += infeasibility; |
2422 | | numberDualInfeasibilities_ ++; |
2423 | | } |
2424 | | kColumn = next_[kColumn]; |
2425 | | } |
2426 | | } |
2427 | | } |
2428 | | } |
2429 | | // and get statistics for column generation |
2430 | | synchronize(model, 4); |
2431 | | infeasibilityWeight_ = -1.0; |
2432 | | } |
2433 | | break; |
2434 | | // Report on infeasibilities of key variables |
2435 | | case 3: { |
2436 | | model->setSumDualInfeasibilities(model->sumDualInfeasibilities() + |
2437 | | sumDualInfeasibilities_); |
2438 | | model->setNumberDualInfeasibilities(model->numberDualInfeasibilities() + |
2439 | | numberDualInfeasibilities_); |
2440 | | model->setSumOfRelaxedDualInfeasibilities(model->sumOfRelaxedDualInfeasibilities() + |
2441 | | sumOfRelaxedDualInfeasibilities_); |
2442 | | } |
2443 | | break; |
2444 | | // modify costs before transposeUpdate for partial pricing |
2445 | | case 4: { |
2446 | | // First compute new costs etc for interesting gubs |
2447 | | int iLook = 0; |
2448 | | int iSet = fromIndex_[0]; |
2449 | | double primalTolerance = model->primalTolerance(); |
2450 | | const double * cost = model->costRegion(); |
2451 | | double * solution = model->solutionRegion(); |
2452 | | double infeasibilityCost = model->infeasibilityCost(); |
2453 | | int numberColumns = model->numberColumns(); |
2454 | | int numberChanged = 0; |
2455 | | int * pivotVariable = model->pivotVariable(); |
2456 | | while (iSet >= 0) { |
2457 | | int key = keyVariable_[iSet]; |
2458 | | double value = 0.0; |
2459 | | // sum over all except key |
2460 | | if ((gubType_ & 8) != 0) { |
2461 | | int iColumn = next_[key]; |
2462 | | // sum all non-key variables |
2463 | | while(iColumn >= 0) { |
2464 | | value += solution[iColumn]; |
2465 | | iColumn = next_[iColumn]; |
2466 | | } |
2467 | | } else { |
2468 | | // bounds exist - sum over all except key |
2469 | | int stop = -(key + 1); |
2470 | | int iColumn = next_[key]; |
2471 | | // sum all non-key variables |
2472 | | while(iColumn != stop) { |
2473 | | if (iColumn < 0) |
2474 | | iColumn = -iColumn - 1; |
2475 | | value += solution[iColumn]; |
2476 | | iColumn = next_[iColumn]; |
2477 | | } |
2478 | | } |
2479 | | double costChange; |
2480 | | double oldCost = changeCost_[iLook]; |
2481 | | if (key < numberColumns) { |
2482 | | assert (getStatus(iSet) != ClpSimplex::basic); |
2483 | | double sol; |
2484 | | if (getStatus(iSet) == ClpSimplex::atUpperBound) |
2485 | | sol = upper_[iSet] - value; |
2486 | | else |
2487 | | sol = lower_[iSet] - value; |
2488 | | solution[key] = sol; |
2489 | | // fix up cost |
2490 | | model->nonLinearCost()->setOne(key, sol); |
2491 | | #ifdef CLP_DEBUG_PRINT |
2492 | | printf("yy Value of key structural %d for set %d is %g - cost %g old cost %g\n", key, iSet, sol, |
2493 | | cost[key], oldCost); |
2494 | | #endif |
2495 | | costChange = cost[key] - oldCost; |
2496 | | } else { |
2497 | | // slack is key |
2498 | | if (value > upper_[iSet] + primalTolerance) { |
2499 | | setAbove(iSet); |
2500 | | } else if (value < lower_[iSet] - primalTolerance) { |
2501 | | setBelow(iSet); |
2502 | | } else { |
2503 | | setFeasible(iSet); |
2504 | | } |
2505 | | // negative as -1.0 for slack |
2506 | | costChange = -weight(iSet) * infeasibilityCost - oldCost; |
2507 | | #ifdef CLP_DEBUG_PRINT |
2508 | | printf("yy Value of key slack for set %d is %g - cost %g old cost %g\n", iSet, value, |
2509 | | weight(iSet)*infeasibilityCost, oldCost); |
2510 | | #endif |
2511 | | } |
2512 | | if (costChange) { |
2513 | | fromIndex_[numberChanged] = iSet; |
2514 | | toIndex_[iSet] = numberChanged; |
2515 | | changeCost_[numberChanged++] = costChange; |
2516 | | } |
2517 | | iSet = fromIndex_[++iLook]; |
2518 | | } |
2519 | | if (numberChanged || possiblePivotKey_ >= 0) { |
2520 | | // first do those in list already |
2521 | | int number = array->getNumElements(); |
2522 | | array->setPacked(); |
2523 | | int i; |
2524 | | double * work = array->denseVector(); |
2525 | | int * which = array->getIndices(); |
2526 | | for (i = 0; i < number; i++) { |
2527 | | int iRow = which[i]; |
2528 | | int iPivot = pivotVariable[iRow]; |
2529 | | if (iPivot < numberColumns) { |
2530 | | int iSet = backward_[iPivot]; |
2531 | | if (iSet >= 0 && toIndex_[iSet] >= 0) { |
2532 | | double newValue = work[i] + changeCost_[toIndex_[iSet]]; |
2533 | | if (!newValue) |
2534 | | newValue = 1.0e-100; |
2535 | | work[i] = newValue; |
2536 | | // mark as done |
2537 | | backward_[iPivot] = -1; |
2538 | | } |
2539 | | } |
2540 | | if (possiblePivotKey_ == iRow) { |
2541 | | double newValue = work[i] - model->dualIn(); |
2542 | | if (!newValue) |
2543 | | newValue = 1.0e-100; |
2544 | | work[i] = newValue; |
2545 | | possiblePivotKey_ = -1; |
2546 | | } |
2547 | | } |
2548 | | // now do rest and clean up |
2549 | | for (i = 0; i < numberChanged; i++) { |
2550 | | int iSet = fromIndex_[i]; |
2551 | | int key = keyVariable_[iSet]; |
2552 | | int iColumn = next_[key]; |
2553 | | double change = changeCost_[i]; |
2554 | | while (iColumn >= 0) { |
2555 | | if (backward_[iColumn] >= 0) { |
2556 | | int iRow = backToPivotRow_[iColumn]; |
2557 | | assert (iRow >= 0); |
2558 | | work[number] = change; |
2559 | | if (possiblePivotKey_ == iRow) { |
2560 | | double newValue = work[number] - model->dualIn(); |
2561 | | if (!newValue) |
2562 | | newValue = 1.0e-100; |
2563 | | work[number] = newValue; |
2564 | | possiblePivotKey_ = -1; |
2565 | | } |
2566 | | which[number++] = iRow; |
2567 | | } else { |
2568 | | // reset |
2569 | | backward_[iColumn] = iSet; |
2570 | | } |
2571 | | iColumn = next_[iColumn]; |
2572 | | } |
2573 | | toIndex_[iSet] = -1; |
2574 | | } |
2575 | | if (possiblePivotKey_ >= 0) { |
2576 | | work[number] = -model->dualIn(); |
2577 | | which[number++] = possiblePivotKey_; |
2578 | | possiblePivotKey_ = -1; |
2579 | | } |
2580 | | fromIndex_[0] = -1; |
2581 | | array->setNumElements(number); |
2582 | | } |
2583 | | } |
2584 | | break; |
2585 | | } |
| 2516 | } |
| 2517 | } |
| 2518 | // now do rest and clean up |
| 2519 | for (i = 0; i < numberChanged; i++) { |
| 2520 | int iSet = fromIndex_[i]; |
| 2521 | int key = keyVariable_[iSet]; |
| 2522 | int iColumn = next_[key]; |
| 2523 | double change = changeCost_[i]; |
| 2524 | while (iColumn >= 0) { |
| 2525 | if (backward_[iColumn] >= 0) { |
| 2526 | int iRow = backToPivotRow_[iColumn]; |
| 2527 | assert(iRow >= 0); |
| 2528 | work[number] = change; |
| 2529 | if (possiblePivotKey_ == iRow) { |
| 2530 | double newValue = work[number] - model->dualIn(); |
| 2531 | if (!newValue) |
| 2532 | newValue = 1.0e-100; |
| 2533 | work[number] = newValue; |
| 2534 | possiblePivotKey_ = -1; |
| 2535 | } |
| 2536 | which[number++] = iRow; |
| 2537 | } else { |
| 2538 | // reset |
| 2539 | backward_[iColumn] = iSet; |
| 2540 | } |
| 2541 | iColumn = next_[iColumn]; |
| 2542 | } |
| 2543 | toIndex_[iSet] = -1; |
| 2544 | } |
| 2545 | if (possiblePivotKey_ >= 0) { |
| 2546 | work[number] = -model->dualIn(); |
| 2547 | which[number++] = possiblePivotKey_; |
| 2548 | possiblePivotKey_ = -1; |
| 2549 | } |
| 2550 | fromIndex_[0] = -1; |
| 2551 | array->setNumElements(number); |
| 2552 | } |
| 2553 | } break; |
| 2554 | } |
2601 | | int returnCode = 0; |
2602 | | int numberColumns = model->numberColumns(); |
2603 | | switch (mode) { |
2604 | | // Fill in pivotVariable but not for key variables |
2605 | | case 0: { |
2606 | | if (!next_ ) { |
2607 | | // do ordering |
2608 | | assert (!rhsOffset_); |
2609 | | // create and do gub crash |
2610 | | useEffectiveRhs(model, false); |
2611 | | } |
2612 | | int i; |
2613 | | int numberBasic = number; |
2614 | | // Use different array so can build from true pivotVariable_ |
2615 | | //int * pivotVariable = model->pivotVariable(); |
2616 | | int * pivotVariable = model->rowArray(0)->getIndices(); |
2617 | | for (i = 0; i < numberColumns; i++) { |
2618 | | if (model->getColumnStatus(i) == ClpSimplex::basic) { |
2619 | | int iSet = backward_[i]; |
2620 | | if (iSet < 0 || i != keyVariable_[iSet]) |
2621 | | pivotVariable[numberBasic++] = i; |
2622 | | } |
2623 | | } |
2624 | | number = numberBasic; |
2625 | | if (model->numberIterations()) |
2626 | | assert (number == model->numberRows()); |
2627 | | } |
2628 | | break; |
2629 | | // Make all key variables basic |
2630 | | case 1: { |
2631 | | int i; |
2632 | | for (i = 0; i < numberSets_; i++) { |
2633 | | int iColumn = keyVariable_[i]; |
2634 | | if (iColumn < numberColumns) |
2635 | | model->setColumnStatus(iColumn, ClpSimplex::basic); |
2636 | | } |
2637 | | } |
2638 | | break; |
2639 | | // Do initial extra rows + maximum basic |
2640 | | case 2: { |
2641 | | returnCode = getNumRows() + 1; |
2642 | | number = model->numberRows() + numberSets_; |
2643 | | } |
2644 | | break; |
2645 | | // Before normal replaceColumn |
2646 | | case 3: { |
2647 | | int sequenceIn = model->sequenceIn(); |
2648 | | int sequenceOut = model->sequenceOut(); |
2649 | | int numberColumns = model->numberColumns(); |
2650 | | int numberRows = model->numberRows(); |
2651 | | int pivotRow = model->pivotRow(); |
2652 | | if (gubSlackIn_ >= 0) |
2653 | | assert (sequenceIn > numberRows + numberColumns); |
2654 | | if (sequenceIn == sequenceOut) |
2655 | | return -1; |
2656 | | int iSetIn = -1; |
2657 | | int iSetOut = -1; |
2658 | | if (sequenceOut < numberColumns) { |
2659 | | iSetOut = backward_[sequenceOut]; |
2660 | | } else if (sequenceOut >= numberRows + numberColumns) { |
2661 | | assert (pivotRow >= numberRows); |
2662 | | int iExtra = pivotRow - numberRows; |
2663 | | assert (iExtra >= 0); |
2664 | | if (iSetOut < 0) |
2665 | | iSetOut = fromIndex_[iExtra]; |
2666 | | else |
2667 | | assert(iSetOut == fromIndex_[iExtra]); |
2668 | | } |
2669 | | if (sequenceIn < numberColumns) { |
2670 | | iSetIn = backward_[sequenceIn]; |
2671 | | } else if (gubSlackIn_ >= 0) { |
2672 | | iSetIn = gubSlackIn_; |
2673 | | } |
2674 | | possiblePivotKey_ = -1; |
2675 | | number = 0; // say do ordinary |
2676 | | int * pivotVariable = model->pivotVariable(); |
2677 | | if (pivotRow >= numberRows) { |
2678 | | int iExtra = pivotRow - numberRows; |
2679 | | //const int * length = matrix_->getVectorLengths(); |
| 2568 | int returnCode = 0; |
| 2569 | int numberColumns = model->numberColumns(); |
| 2570 | switch (mode) { |
| 2571 | // Fill in pivotVariable but not for key variables |
| 2572 | case 0: { |
| 2573 | if (!next_) { |
| 2574 | // do ordering |
| 2575 | assert(!rhsOffset_); |
| 2576 | // create and do gub crash |
| 2577 | useEffectiveRhs(model, false); |
| 2578 | } |
| 2579 | int i; |
| 2580 | int numberBasic = number; |
| 2581 | // Use different array so can build from true pivotVariable_ |
| 2582 | //int * pivotVariable = model->pivotVariable(); |
| 2583 | int *pivotVariable = model->rowArray(0)->getIndices(); |
| 2584 | for (i = 0; i < numberColumns; i++) { |
| 2585 | if (model->getColumnStatus(i) == ClpSimplex::basic) { |
| 2586 | int iSet = backward_[i]; |
| 2587 | if (iSet < 0 || i != keyVariable_[iSet]) |
| 2588 | pivotVariable[numberBasic++] = i; |
| 2589 | } |
| 2590 | } |
| 2591 | number = numberBasic; |
| 2592 | if (model->numberIterations()) |
| 2593 | assert(number == model->numberRows()); |
| 2594 | } break; |
| 2595 | // Make all key variables basic |
| 2596 | case 1: { |
| 2597 | int i; |
| 2598 | for (i = 0; i < numberSets_; i++) { |
| 2599 | int iColumn = keyVariable_[i]; |
| 2600 | if (iColumn < numberColumns) |
| 2601 | model->setColumnStatus(iColumn, ClpSimplex::basic); |
| 2602 | } |
| 2603 | } break; |
| 2604 | // Do initial extra rows + maximum basic |
| 2605 | case 2: { |
| 2606 | returnCode = getNumRows() + 1; |
| 2607 | number = model->numberRows() + numberSets_; |
| 2608 | } break; |
| 2609 | // Before normal replaceColumn |
| 2610 | case 3: { |
| 2611 | int sequenceIn = model->sequenceIn(); |
| 2612 | int sequenceOut = model->sequenceOut(); |
| 2613 | int numberColumns = model->numberColumns(); |
| 2614 | int numberRows = model->numberRows(); |
| 2615 | int pivotRow = model->pivotRow(); |
| 2616 | if (gubSlackIn_ >= 0) |
| 2617 | assert(sequenceIn > numberRows + numberColumns); |
| 2618 | if (sequenceIn == sequenceOut) |
| 2619 | return -1; |
| 2620 | int iSetIn = -1; |
| 2621 | int iSetOut = -1; |
| 2622 | if (sequenceOut < numberColumns) { |
| 2623 | iSetOut = backward_[sequenceOut]; |
| 2624 | } else if (sequenceOut >= numberRows + numberColumns) { |
| 2625 | assert(pivotRow >= numberRows); |
| 2626 | int iExtra = pivotRow - numberRows; |
| 2627 | assert(iExtra >= 0); |
| 2628 | if (iSetOut < 0) |
| 2629 | iSetOut = fromIndex_[iExtra]; |
| 2630 | else |
| 2631 | assert(iSetOut == fromIndex_[iExtra]); |
| 2632 | } |
| 2633 | if (sequenceIn < numberColumns) { |
| 2634 | iSetIn = backward_[sequenceIn]; |
| 2635 | } else if (gubSlackIn_ >= 0) { |
| 2636 | iSetIn = gubSlackIn_; |
| 2637 | } |
| 2638 | possiblePivotKey_ = -1; |
| 2639 | number = 0; // say do ordinary |
| 2640 | int *pivotVariable = model->pivotVariable(); |
| 2641 | if (pivotRow >= numberRows) { |
| 2642 | int iExtra = pivotRow - numberRows; |
| 2643 | //const int * length = matrix_->getVectorLengths(); |