73 | | otherFlags_ = rhs.otherFlags_; |
74 | | #endif |
75 | | columnOrdered_ = rhs.columnOrdered_; |
76 | | if (numberColumns_) { |
77 | | CoinBigIndex numberElements = rhs.startPositive_[numberColumns_]; |
78 | | indices_ = new int [ numberElements]; |
79 | | CoinMemcpyN(rhs.indices_, numberElements, indices_); |
80 | | startPositive_ = new CoinBigIndex [ numberColumns_+1]; |
81 | | CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_); |
82 | | startNegative_ = new CoinBigIndex [ numberColumns_]; |
83 | | CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_); |
84 | | } |
85 | | int numberRows = getNumRows(); |
86 | | if (rhs.rhsOffset_ && numberRows) { |
87 | | rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows); |
88 | | } else { |
89 | | rhsOffset_ = NULL; |
90 | | } |
| 72 | otherFlags_ = rhs.otherFlags_; |
| 73 | #endif |
| 74 | columnOrdered_ = rhs.columnOrdered_; |
| 75 | if (numberColumns_) { |
| 76 | CoinBigIndex numberElements = rhs.startPositive_[numberColumns_]; |
| 77 | indices_ = new int[numberElements]; |
| 78 | CoinMemcpyN(rhs.indices_, numberElements, indices_); |
| 79 | startPositive_ = new CoinBigIndex[numberColumns_ + 1]; |
| 80 | CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_); |
| 81 | startNegative_ = new CoinBigIndex[numberColumns_]; |
| 82 | CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_); |
| 83 | } |
| 84 | int numberRows = getNumRows(); |
| 85 | if (rhs.rhsOffset_ && numberRows) { |
| 86 | rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows); |
| 87 | } else { |
| 88 | rhsOffset_ = NULL; |
| 89 | } |
117 | | ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const CoinPackedMatrix & rhs) |
118 | | : ClpMatrixBase() |
119 | | { |
120 | | setType(12); |
121 | | matrix_ = NULL; |
122 | | startPositive_ = NULL; |
123 | | startNegative_ = NULL; |
124 | | lengths_ = NULL; |
125 | | indices_ = NULL; |
126 | | int iColumn; |
127 | | assert (rhs.isColOrdered()); |
128 | | // get matrix data pointers |
129 | | const int * row = rhs.getIndices(); |
130 | | const CoinBigIndex * columnStart = rhs.getVectorStarts(); |
131 | | const int * columnLength = rhs.getVectorLengths(); |
132 | | const double * elementByColumn = rhs.getElements(); |
133 | | numberColumns_ = rhs.getNumCols(); |
134 | | numberRows_ = -1; |
135 | | indices_ = new int[rhs.getNumElements()]; |
136 | | startPositive_ = new CoinBigIndex [numberColumns_+1]; |
137 | | startNegative_ = new CoinBigIndex [numberColumns_]; |
138 | | int * temp = new int [rhs.getNumRows()]; |
139 | | CoinBigIndex j = 0; |
140 | | CoinBigIndex numberGoodP = 0; |
141 | | CoinBigIndex numberGoodM = 0; |
142 | | CoinBigIndex numberBad = 0; |
143 | | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
144 | | CoinBigIndex k; |
145 | | int iNeg = 0; |
146 | | startPositive_[iColumn] = j; |
147 | | for (k = columnStart[iColumn]; k < columnStart[iColumn] + columnLength[iColumn]; |
148 | | k++) { |
149 | | int iRow; |
150 | | if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) { |
151 | | iRow = row[k]; |
152 | | numberRows_ = CoinMax(numberRows_, iRow); |
153 | | indices_[j++] = iRow; |
154 | | numberGoodP++; |
155 | | } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) { |
156 | | iRow = row[k]; |
157 | | numberRows_ = CoinMax(numberRows_, iRow); |
158 | | temp[iNeg++] = iRow; |
159 | | numberGoodM++; |
160 | | } else { |
161 | | numberBad++; |
162 | | } |
163 | | } |
164 | | // move negative |
165 | | startNegative_[iColumn] = j; |
166 | | for (k = 0; k < iNeg; k++) { |
167 | | indices_[j++] = temp[k]; |
168 | | } |
169 | | } |
170 | | startPositive_[numberColumns_] = j; |
171 | | delete [] temp; |
172 | | if (numberBad) { |
173 | | delete [] indices_; |
174 | | indices_ = NULL; |
175 | | numberRows_ = 0; |
176 | | numberColumns_ = 0; |
177 | | delete [] startPositive_; |
178 | | delete [] startNegative_; |
179 | | // Put in statistics |
180 | | startPositive_ = new CoinBigIndex [3]; |
181 | | startPositive_[0] = numberGoodP; |
182 | | startPositive_[1] = numberGoodM; |
183 | | startPositive_[2] = numberBad; |
184 | | startNegative_ = NULL; |
185 | | } else { |
186 | | numberRows_ ++; // correct |
187 | | // but number should be same as rhs |
188 | | assert (numberRows_ <= rhs.getNumRows()); |
189 | | numberRows_ = rhs.getNumRows(); |
190 | | columnOrdered_ = true; |
191 | | } |
192 | | // Check valid |
193 | | if (!numberBad) |
194 | | checkValid(false); |
| 116 | ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix(const CoinPackedMatrix &rhs) |
| 117 | : ClpMatrixBase() |
| 118 | { |
| 119 | setType(12); |
| 120 | matrix_ = NULL; |
| 121 | startPositive_ = NULL; |
| 122 | startNegative_ = NULL; |
| 123 | lengths_ = NULL; |
| 124 | indices_ = NULL; |
| 125 | int iColumn; |
| 126 | assert(rhs.isColOrdered()); |
| 127 | // get matrix data pointers |
| 128 | const int *row = rhs.getIndices(); |
| 129 | const CoinBigIndex *columnStart = rhs.getVectorStarts(); |
| 130 | const int *columnLength = rhs.getVectorLengths(); |
| 131 | const double *elementByColumn = rhs.getElements(); |
| 132 | numberColumns_ = rhs.getNumCols(); |
| 133 | numberRows_ = -1; |
| 134 | indices_ = new int[rhs.getNumElements()]; |
| 135 | startPositive_ = new CoinBigIndex[numberColumns_ + 1]; |
| 136 | startNegative_ = new CoinBigIndex[numberColumns_]; |
| 137 | int *temp = new int[rhs.getNumRows()]; |
| 138 | CoinBigIndex j = 0; |
| 139 | CoinBigIndex numberGoodP = 0; |
| 140 | CoinBigIndex numberGoodM = 0; |
| 141 | CoinBigIndex numberBad = 0; |
| 142 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
| 143 | CoinBigIndex k; |
| 144 | int iNeg = 0; |
| 145 | startPositive_[iColumn] = j; |
| 146 | for (k = columnStart[iColumn]; k < columnStart[iColumn] + columnLength[iColumn]; |
| 147 | k++) { |
| 148 | int iRow; |
| 149 | if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) { |
| 150 | iRow = row[k]; |
| 151 | numberRows_ = CoinMax(numberRows_, iRow); |
| 152 | indices_[j++] = iRow; |
| 153 | numberGoodP++; |
| 154 | } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) { |
| 155 | iRow = row[k]; |
| 156 | numberRows_ = CoinMax(numberRows_, iRow); |
| 157 | temp[iNeg++] = iRow; |
| 158 | numberGoodM++; |
| 159 | } else { |
| 160 | numberBad++; |
| 161 | } |
| 162 | } |
| 163 | // move negative |
| 164 | startNegative_[iColumn] = j; |
| 165 | for (k = 0; k < iNeg; k++) { |
| 166 | indices_[j++] = temp[k]; |
| 167 | } |
| 168 | } |
| 169 | startPositive_[numberColumns_] = j; |
| 170 | delete[] temp; |
| 171 | if (numberBad) { |
| 172 | delete[] indices_; |
| 173 | indices_ = NULL; |
| 174 | numberRows_ = 0; |
| 175 | numberColumns_ = 0; |
| 176 | delete[] startPositive_; |
| 177 | delete[] startNegative_; |
| 178 | // Put in statistics |
| 179 | startPositive_ = new CoinBigIndex[3]; |
| 180 | startPositive_[0] = numberGoodP; |
| 181 | startPositive_[1] = numberGoodM; |
| 182 | startPositive_[2] = numberBad; |
| 183 | startNegative_ = NULL; |
| 184 | } else { |
| 185 | numberRows_++; // correct |
| 186 | // but number should be same as rhs |
| 187 | assert(numberRows_ <= rhs.getNumRows()); |
| 188 | numberRows_ = rhs.getNumRows(); |
| 189 | columnOrdered_ = true; |
| 190 | } |
| 191 | // Check valid |
| 192 | if (!numberBad) |
| 193 | checkValid(false); |
230 | | otherFlags_ = rhs.otherFlags_; |
231 | | #endif |
232 | | if (numberColumns_) { |
233 | | CoinBigIndex numberElements = rhs.startPositive_[numberColumns_]; |
234 | | indices_ = new int [ numberElements]; |
235 | | CoinMemcpyN(rhs.indices_, numberElements, indices_); |
236 | | startPositive_ = new CoinBigIndex [ numberColumns_+1]; |
237 | | CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_); |
238 | | startNegative_ = new CoinBigIndex [ numberColumns_]; |
239 | | CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_); |
240 | | } |
241 | | } |
242 | | return *this; |
| 229 | otherFlags_ = rhs.otherFlags_; |
| 230 | #endif |
| 231 | if (numberColumns_) { |
| 232 | CoinBigIndex numberElements = rhs.startPositive_[numberColumns_]; |
| 233 | indices_ = new int[numberElements]; |
| 234 | CoinMemcpyN(rhs.indices_, numberElements, indices_); |
| 235 | startPositive_ = new CoinBigIndex[numberColumns_ + 1]; |
| 236 | CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_); |
| 237 | startNegative_ = new CoinBigIndex[numberColumns_]; |
| 238 | CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_); |
| 239 | } |
| 240 | } |
| 241 | return *this; |
289 | | int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_; |
290 | | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
291 | | int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_; |
292 | | int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_; |
293 | | // Also swap incoming if not column ordered |
294 | | if (!columnOrdered_) { |
295 | | int temp1 = numberRows; |
296 | | numberRows = numberColumns; |
297 | | numberColumns = temp1; |
298 | | const int * temp2; |
299 | | temp2 = whichRow; |
300 | | whichRow = whichColumn; |
301 | | whichColumn = temp2; |
| 288 | int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_; |
| 289 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
| 290 | int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_; |
| 291 | int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_; |
| 292 | // Also swap incoming if not column ordered |
| 293 | if (!columnOrdered_) { |
| 294 | int temp1 = numberRows; |
| 295 | numberRows = numberColumns; |
| 296 | numberColumns = temp1; |
| 297 | const int *temp2; |
| 298 | temp2 = whichRow; |
| 299 | whichRow = whichColumn; |
| 300 | whichColumn = temp2; |
| 301 | } |
| 302 | // Throw exception if rhs empty |
| 303 | if (numberMajor1 <= 0 || numberMinor1 <= 0) |
| 304 | throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix"); |
| 305 | // Array to say if an old row is in new copy |
| 306 | int *newRow = new int[numberMinor1]; |
| 307 | int iRow; |
| 308 | for (iRow = 0; iRow < numberMinor1; iRow++) |
| 309 | newRow[iRow] = -1; |
| 310 | // and array for duplicating rows |
| 311 | int *duplicateRow = new int[numberMinor]; |
| 312 | int numberBad = 0; |
| 313 | for (iRow = 0; iRow < numberMinor; iRow++) { |
| 314 | duplicateRow[iRow] = -1; |
| 315 | int kRow = whichRow[iRow]; |
| 316 | if (kRow >= 0 && kRow < numberMinor1) { |
| 317 | if (newRow[kRow] < 0) { |
| 318 | // first time |
| 319 | newRow[kRow] = iRow; |
| 320 | } else { |
| 321 | // duplicate |
| 322 | int lastRow = newRow[kRow]; |
| 323 | newRow[kRow] = iRow; |
| 324 | duplicateRow[iRow] = lastRow; |
| 325 | } |
| 326 | } else { |
| 327 | // bad row |
| 328 | numberBad++; |
| 329 | } |
| 330 | } |
| 331 | |
| 332 | if (numberBad) |
| 333 | throw CoinError("bad minor entries", |
| 334 | "subset constructor", "ClpPlusMinusOneMatrix"); |
| 335 | // now get size and check columns |
| 336 | CoinBigIndex size = 0; |
| 337 | int iColumn; |
| 338 | numberBad = 0; |
| 339 | for (iColumn = 0; iColumn < numberMajor; iColumn++) { |
| 340 | int kColumn = whichColumn[iColumn]; |
| 341 | if (kColumn >= 0 && kColumn < numberMajor1) { |
| 342 | CoinBigIndex i; |
| 343 | for (i = startPositive1[kColumn]; i < startPositive1[kColumn + 1]; i++) { |
| 344 | int kRow = index1[i]; |
| 345 | kRow = newRow[kRow]; |
| 346 | while (kRow >= 0) { |
| 347 | size++; |
| 348 | kRow = duplicateRow[kRow]; |
303 | | // Throw exception if rhs empty |
304 | | if (numberMajor1 <= 0 || numberMinor1 <= 0) |
305 | | throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix"); |
306 | | // Array to say if an old row is in new copy |
307 | | int * newRow = new int [numberMinor1]; |
308 | | int iRow; |
309 | | for (iRow = 0; iRow < numberMinor1; iRow++) |
310 | | newRow[iRow] = -1; |
311 | | // and array for duplicating rows |
312 | | int * duplicateRow = new int [numberMinor]; |
313 | | int numberBad = 0; |
314 | | for (iRow = 0; iRow < numberMinor; iRow++) { |
315 | | duplicateRow[iRow] = -1; |
316 | | int kRow = whichRow[iRow]; |
317 | | if (kRow >= 0 && kRow < numberMinor1) { |
318 | | if (newRow[kRow] < 0) { |
319 | | // first time |
320 | | newRow[kRow] = iRow; |
321 | | } else { |
322 | | // duplicate |
323 | | int lastRow = newRow[kRow]; |
324 | | newRow[kRow] = iRow; |
325 | | duplicateRow[iRow] = lastRow; |
326 | | } |
327 | | } else { |
328 | | // bad row |
329 | | numberBad++; |
330 | | } |
331 | | } |
332 | | |
333 | | if (numberBad) |
334 | | throw CoinError("bad minor entries", |
335 | | "subset constructor", "ClpPlusMinusOneMatrix"); |
336 | | // now get size and check columns |
337 | | CoinBigIndex size = 0; |
338 | | int iColumn; |
339 | | numberBad = 0; |
340 | | for (iColumn = 0; iColumn < numberMajor; iColumn++) { |
341 | | int kColumn = whichColumn[iColumn]; |
342 | | if (kColumn >= 0 && kColumn < numberMajor1) { |
343 | | CoinBigIndex i; |
344 | | for (i = startPositive1[kColumn]; i < startPositive1[kColumn+1]; i++) { |
345 | | int kRow = index1[i]; |
346 | | kRow = newRow[kRow]; |
347 | | while (kRow >= 0) { |
348 | | size++; |
349 | | kRow = duplicateRow[kRow]; |
350 | | } |
351 | | } |
352 | | } else { |
353 | | // bad column |
354 | | numberBad++; |
355 | | printf("%d %d %d %d\n", iColumn, numberMajor, numberMajor1, kColumn); |
356 | | } |
357 | | } |
358 | | if (numberBad) |
359 | | throw CoinError("bad major entries", |
360 | | "subset constructor", "ClpPlusMinusOneMatrix"); |
361 | | // now create arrays |
362 | | startPositive_ = new CoinBigIndex [numberMajor+1]; |
363 | | startNegative_ = new CoinBigIndex [numberMajor]; |
364 | | indices_ = new int[size]; |
365 | | // and fill them |
366 | | size = 0; |
367 | | startPositive_[0] = 0; |
368 | | CoinBigIndex * startNegative1 = rhs.startNegative_; |
369 | | for (iColumn = 0; iColumn < numberMajor; iColumn++) { |
370 | | int kColumn = whichColumn[iColumn]; |
371 | | CoinBigIndex i; |
372 | | for (i = startPositive1[kColumn]; i < startNegative1[kColumn]; i++) { |
373 | | int kRow = index1[i]; |
374 | | kRow = newRow[kRow]; |
375 | | while (kRow >= 0) { |
376 | | indices_[size++] = kRow; |
377 | | kRow = duplicateRow[kRow]; |
378 | | } |
379 | | } |
380 | | startNegative_[iColumn] = size; |
381 | | for (; i < startPositive1[kColumn+1]; i++) { |
382 | | int kRow = index1[i]; |
383 | | kRow = newRow[kRow]; |
384 | | while (kRow >= 0) { |
385 | | indices_[size++] = kRow; |
386 | | kRow = duplicateRow[kRow]; |
387 | | } |
388 | | } |
389 | | startPositive_[iColumn+1] = size; |
390 | | } |
391 | | delete [] newRow; |
392 | | delete [] duplicateRow; |
393 | | } |
394 | | // Check valid |
395 | | checkValid(false); |
396 | | } |
397 | | |
| 350 | } |
| 351 | } else { |
| 352 | // bad column |
| 353 | numberBad++; |
| 354 | printf("%d %d %d %d\n", iColumn, numberMajor, numberMajor1, kColumn); |
| 355 | } |
| 356 | } |
| 357 | if (numberBad) |
| 358 | throw CoinError("bad major entries", |
| 359 | "subset constructor", "ClpPlusMinusOneMatrix"); |
| 360 | // now create arrays |
| 361 | startPositive_ = new CoinBigIndex[numberMajor + 1]; |
| 362 | startNegative_ = new CoinBigIndex[numberMajor]; |
| 363 | indices_ = new int[size]; |
| 364 | // and fill them |
| 365 | size = 0; |
| 366 | startPositive_[0] = 0; |
| 367 | CoinBigIndex *startNegative1 = rhs.startNegative_; |
| 368 | for (iColumn = 0; iColumn < numberMajor; iColumn++) { |
| 369 | int kColumn = whichColumn[iColumn]; |
| 370 | CoinBigIndex i; |
| 371 | for (i = startPositive1[kColumn]; i < startNegative1[kColumn]; i++) { |
| 372 | int kRow = index1[i]; |
| 373 | kRow = newRow[kRow]; |
| 374 | while (kRow >= 0) { |
| 375 | indices_[size++] = kRow; |
| 376 | kRow = duplicateRow[kRow]; |
| 377 | } |
| 378 | } |
| 379 | startNegative_[iColumn] = size; |
| 380 | for (; i < startPositive1[kColumn + 1]; i++) { |
| 381 | int kRow = index1[i]; |
| 382 | kRow = newRow[kRow]; |
| 383 | while (kRow >= 0) { |
| 384 | indices_[size++] = kRow; |
| 385 | kRow = duplicateRow[kRow]; |
| 386 | } |
| 387 | } |
| 388 | startPositive_[iColumn + 1] = size; |
| 389 | } |
| 390 | delete[] newRow; |
| 391 | delete[] duplicateRow; |
| 392 | } |
| 393 | // Check valid |
| 394 | checkValid(false); |
| 395 | } |
403 | | int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_; |
404 | | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
405 | | // count number in each row/column |
406 | | CoinBigIndex * tempP = new CoinBigIndex [numberMinor]; |
407 | | CoinBigIndex * tempN = new CoinBigIndex [numberMinor]; |
408 | | memset(tempP, 0, numberMinor * sizeof(CoinBigIndex)); |
409 | | memset(tempN, 0, numberMinor * sizeof(CoinBigIndex)); |
410 | | CoinBigIndex j = 0; |
411 | | int i; |
412 | | for (i = 0; i < numberMajor; i++) { |
413 | | for (; j < startNegative_[i]; j++) { |
414 | | int iRow = indices_[j]; |
415 | | tempP[iRow]++; |
416 | | } |
417 | | for (; j < startPositive_[i+1]; j++) { |
418 | | int iRow = indices_[j]; |
419 | | tempN[iRow]++; |
420 | | } |
421 | | } |
422 | | int * newIndices = new int [startPositive_[numberMajor]]; |
423 | | CoinBigIndex * newP = new CoinBigIndex [numberMinor+1]; |
424 | | CoinBigIndex * newN = new CoinBigIndex[numberMinor]; |
425 | | int iRow; |
426 | | j = 0; |
427 | | // do starts |
428 | | for (iRow = 0; iRow < numberMinor; iRow++) { |
429 | | newP[iRow] = j; |
430 | | j += tempP[iRow]; |
431 | | tempP[iRow] = newP[iRow]; |
432 | | newN[iRow] = j; |
433 | | j += tempN[iRow]; |
434 | | tempN[iRow] = newN[iRow]; |
435 | | } |
436 | | newP[numberMinor] = j; |
437 | | j = 0; |
438 | | for (i = 0; i < numberMajor; i++) { |
439 | | for (; j < startNegative_[i]; j++) { |
440 | | int iRow = indices_[j]; |
441 | | CoinBigIndex put = tempP[iRow]; |
442 | | newIndices[put++] = i; |
443 | | tempP[iRow] = put; |
444 | | } |
445 | | for (; j < startPositive_[i+1]; j++) { |
446 | | int iRow = indices_[j]; |
447 | | CoinBigIndex put = tempN[iRow]; |
448 | | newIndices[put++] = i; |
449 | | tempN[iRow] = put; |
450 | | } |
451 | | } |
452 | | delete [] tempP; |
453 | | delete [] tempN; |
454 | | ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix(); |
455 | | newCopy->passInCopy(numberMinor, numberMajor, |
456 | | !columnOrdered_, newIndices, newP, newN); |
457 | | return newCopy; |
| 401 | int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_; |
| 402 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
| 403 | // count number in each row/column |
| 404 | CoinBigIndex *tempP = new CoinBigIndex[numberMinor]; |
| 405 | CoinBigIndex *tempN = new CoinBigIndex[numberMinor]; |
| 406 | memset(tempP, 0, numberMinor * sizeof(CoinBigIndex)); |
| 407 | memset(tempN, 0, numberMinor * sizeof(CoinBigIndex)); |
| 408 | CoinBigIndex j = 0; |
| 409 | int i; |
| 410 | for (i = 0; i < numberMajor; i++) { |
| 411 | for (; j < startNegative_[i]; j++) { |
| 412 | int iRow = indices_[j]; |
| 413 | tempP[iRow]++; |
| 414 | } |
| 415 | for (; j < startPositive_[i + 1]; j++) { |
| 416 | int iRow = indices_[j]; |
| 417 | tempN[iRow]++; |
| 418 | } |
| 419 | } |
| 420 | int *newIndices = new int[startPositive_[numberMajor]]; |
| 421 | CoinBigIndex *newP = new CoinBigIndex[numberMinor + 1]; |
| 422 | CoinBigIndex *newN = new CoinBigIndex[numberMinor]; |
| 423 | int iRow; |
| 424 | j = 0; |
| 425 | // do starts |
| 426 | for (iRow = 0; iRow < numberMinor; iRow++) { |
| 427 | newP[iRow] = j; |
| 428 | j += tempP[iRow]; |
| 429 | tempP[iRow] = newP[iRow]; |
| 430 | newN[iRow] = j; |
| 431 | j += tempN[iRow]; |
| 432 | tempN[iRow] = newN[iRow]; |
| 433 | } |
| 434 | newP[numberMinor] = j; |
| 435 | j = 0; |
| 436 | for (i = 0; i < numberMajor; i++) { |
| 437 | for (; j < startNegative_[i]; j++) { |
| 438 | int iRow = indices_[j]; |
| 439 | CoinBigIndex put = tempP[iRow]; |
| 440 | newIndices[put++] = i; |
| 441 | tempP[iRow] = put; |
| 442 | } |
| 443 | for (; j < startPositive_[i + 1]; j++) { |
| 444 | int iRow = indices_[j]; |
| 445 | CoinBigIndex put = tempN[iRow]; |
| 446 | newIndices[put++] = i; |
| 447 | tempN[iRow] = put; |
| 448 | } |
| 449 | } |
| 450 | delete[] tempP; |
| 451 | delete[] tempN; |
| 452 | ClpPlusMinusOneMatrix *newCopy = new ClpPlusMinusOneMatrix(); |
| 453 | newCopy->passInCopy(numberMinor, numberMajor, |
| 454 | !columnOrdered_, newIndices, newP, newN); |
| 455 | return newCopy; |
525 | | } else { |
526 | | // plus one |
527 | | oneit(1); |
528 | | for (i = 0; i < numberMajor; i++) { |
529 | | double value = 0.0; |
530 | | for (; j < startPositive_[i+1]; j++) { |
531 | | int iRow = indices_[j]; |
532 | | value += x[iRow]; |
533 | | } |
534 | | y[i] += scalar * value; |
535 | | } |
536 | | } |
537 | | #endif |
538 | | } |
539 | | void |
540 | | ClpPlusMinusOneMatrix::times(double scalar, |
541 | | const double * COIN_RESTRICT x, double * COIN_RESTRICT y, |
542 | | const double * /*rowScale*/, |
543 | | const double * /*columnScale*/) const |
544 | | { |
545 | | // we know it is not scaled |
546 | | times(scalar, x, y); |
547 | | } |
548 | | void |
549 | | ClpPlusMinusOneMatrix::transposeTimes( double scalar, |
550 | | const double * COIN_RESTRICT x, double * COIN_RESTRICT y, |
551 | | const double * /*rowScale*/, |
552 | | const double * /*columnScale*/, |
553 | | double * /*spare*/) const |
554 | | { |
555 | | // we know it is not scaled |
556 | | transposeTimes(scalar, x, y); |
| 521 | } else { |
| 522 | // plus one |
| 523 | oneit(1); |
| 524 | for (i = 0; i < numberMajor; i++) { |
| 525 | double value = 0.0; |
| 526 | for (; j < startPositive_[i + 1]; j++) { |
| 527 | int iRow = indices_[j]; |
| 528 | value += x[iRow]; |
| 529 | } |
| 530 | y[i] += scalar * value; |
| 531 | } |
| 532 | } |
| 533 | #endif |
| 534 | } |
| 535 | void ClpPlusMinusOneMatrix::times(double scalar, |
| 536 | const double *COIN_RESTRICT x, double *COIN_RESTRICT y, |
| 537 | const double * /*rowScale*/, |
| 538 | const double * /*columnScale*/) const |
| 539 | { |
| 540 | // we know it is not scaled |
| 541 | times(scalar, x, y); |
| 542 | } |
| 543 | void ClpPlusMinusOneMatrix::transposeTimes(double scalar, |
| 544 | const double *COIN_RESTRICT x, double *COIN_RESTRICT y, |
| 545 | const double * /*rowScale*/, |
| 546 | const double * /*columnScale*/, |
| 547 | double * /*spare*/) const |
| 548 | { |
| 549 | // we know it is not scaled |
| 550 | transposeTimes(scalar, x, y); |
560 | | void |
561 | | ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex * model, double scalar, |
562 | | const CoinIndexedVector * rowArray, |
563 | | CoinIndexedVector * y, |
564 | | CoinIndexedVector * columnArray) const |
565 | | { |
566 | | // we know it is not scaled |
567 | | columnArray->clear(); |
568 | | double * COIN_RESTRICT pi = rowArray->denseVector(); |
569 | | int numberNonZero = 0; |
570 | | int * COIN_RESTRICT index = columnArray->getIndices(); |
571 | | double * COIN_RESTRICT array = columnArray->denseVector(); |
572 | | int numberInRowArray = rowArray->getNumElements(); |
573 | | // maybe I need one in OsiSimplex |
574 | | double zeroTolerance = model->zeroTolerance(); |
575 | | int numberRows = model->numberRows(); |
576 | | bool packed = rowArray->packedMode(); |
| 554 | void ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex *model, double scalar, |
| 555 | const CoinIndexedVector *rowArray, |
| 556 | CoinIndexedVector *y, |
| 557 | CoinIndexedVector *columnArray) const |
| 558 | { |
| 559 | // we know it is not scaled |
| 560 | columnArray->clear(); |
| 561 | double *COIN_RESTRICT pi = rowArray->denseVector(); |
| 562 | int numberNonZero = 0; |
| 563 | int *COIN_RESTRICT index = columnArray->getIndices(); |
| 564 | double *COIN_RESTRICT array = columnArray->denseVector(); |
| 565 | int numberInRowArray = rowArray->getNumElements(); |
| 566 | // maybe I need one in OsiSimplex |
| 567 | double zeroTolerance = model->zeroTolerance(); |
| 568 | int numberRows = model->numberRows(); |
| 569 | bool packed = rowArray->packedMode(); |
581 | | ClpPlusMinusOneMatrix* rowCopy = |
582 | | static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy()); |
583 | | #endif |
584 | | double factor = 0.3; |
585 | | // We may not want to do by row if there may be cache problems |
586 | | int numberColumns = model->numberColumns(); |
587 | | // It would be nice to find L2 cache size - for moment 512K |
588 | | // Be slightly optimistic |
589 | | if (numberColumns * sizeof(double) > 1000000) { |
590 | | if (numberRows * 10 < numberColumns) |
591 | | factor = 0.1; |
592 | | else if (numberRows * 4 < numberColumns) |
593 | | factor = 0.15; |
594 | | else if (numberRows * 2 < numberColumns) |
595 | | factor = 0.2; |
596 | | } |
597 | | if (numberInRowArray > factor * numberRows || !rowCopy) { |
598 | | assert (!y->getNumElements()); |
599 | | // do by column |
600 | | // Need to expand if packed mode |
601 | | int iColumn; |
602 | | CoinBigIndex j = 0; |
603 | | assert (columnOrdered_); |
604 | | if (packed) { |
605 | | // need to expand pi into y |
606 | | assert(y->capacity() >= numberRows); |
607 | | double * COIN_RESTRICT piOld = pi; |
608 | | pi = y->denseVector(); |
609 | | const int * COIN_RESTRICT whichRow = rowArray->getIndices(); |
610 | | int i; |
611 | | // modify pi so can collapse to one loop |
612 | | for (i = 0; i < numberInRowArray; i++) { |
613 | | int iRow = whichRow[i]; |
614 | | pi[iRow] = scalar * piOld[i]; |
615 | | } |
| 573 | ClpPlusMinusOneMatrix *rowCopy = static_cast< ClpPlusMinusOneMatrix * >(model->rowCopy()); |
| 574 | #endif |
| 575 | double factor = 0.3; |
| 576 | // We may not want to do by row if there may be cache problems |
| 577 | int numberColumns = model->numberColumns(); |
| 578 | // It would be nice to find L2 cache size - for moment 512K |
| 579 | // Be slightly optimistic |
| 580 | if (numberColumns * sizeof(double) > 1000000) { |
| 581 | if (numberRows * 10 < numberColumns) |
| 582 | factor = 0.1; |
| 583 | else if (numberRows * 4 < numberColumns) |
| 584 | factor = 0.15; |
| 585 | else if (numberRows * 2 < numberColumns) |
| 586 | factor = 0.2; |
| 587 | } |
| 588 | if (numberInRowArray > factor * numberRows || !rowCopy) { |
| 589 | assert(!y->getNumElements()); |
| 590 | // do by column |
| 591 | // Need to expand if packed mode |
| 592 | int iColumn; |
| 593 | CoinBigIndex j = 0; |
| 594 | assert(columnOrdered_); |
| 595 | if (packed) { |
| 596 | // need to expand pi into y |
| 597 | assert(y->capacity() >= numberRows); |
| 598 | double *COIN_RESTRICT piOld = pi; |
| 599 | pi = y->denseVector(); |
| 600 | const int *COIN_RESTRICT whichRow = rowArray->getIndices(); |
| 601 | int i; |
| 602 | // modify pi so can collapse to one loop |
| 603 | for (i = 0; i < numberInRowArray; i++) { |
| 604 | int iRow = whichRow[i]; |
| 605 | pi[iRow] = scalar * piOld[i]; |
| 606 | } |
681 | | void |
682 | | ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar, |
683 | | const CoinIndexedVector * rowArray, |
684 | | CoinIndexedVector * y, |
685 | | CoinIndexedVector * columnArray) const |
686 | | { |
687 | | columnArray->clear(); |
688 | | double * COIN_RESTRICT pi = rowArray->denseVector(); |
689 | | int numberNonZero = 0; |
690 | | int * COIN_RESTRICT index = columnArray->getIndices(); |
691 | | double * COIN_RESTRICT array = columnArray->denseVector(); |
692 | | int numberInRowArray = rowArray->getNumElements(); |
693 | | // maybe I need one in OsiSimplex |
694 | | double zeroTolerance = model->zeroTolerance(); |
695 | | const int * COIN_RESTRICT column = indices_; |
696 | | const CoinBigIndex * COIN_RESTRICT startPositive = startPositive_; |
697 | | const CoinBigIndex * COIN_RESTRICT startNegative = startNegative_; |
698 | | const int * COIN_RESTRICT whichRow = rowArray->getIndices(); |
699 | | bool packed = rowArray->packedMode(); |
700 | | if (numberInRowArray > 2) { |
701 | | // do by rows |
702 | | int iRow; |
703 | | double * COIN_RESTRICT markVector = y->denseVector(); // probably empty .. but |
704 | | int numberOriginal = 0; |
705 | | int i; |
706 | | if (packed) { |
707 | | CoinBigIndex numberCovered=0; |
708 | | int numberColumns = getNumCols(); |
709 | | bool sparse=true; |
710 | | int target=1*numberColumns; |
711 | | for (int i = 0; i < numberInRowArray; i++) { |
712 | | int iRow = whichRow[i]; |
713 | | numberCovered += startPositive[iRow+1] - startPositive[iRow]; |
714 | | if (numberCovered>target) { |
715 | | sparse=false; |
716 | | break; |
717 | | } |
718 | | } |
719 | | numberNonZero = 0; |
720 | | if (sparse) { |
721 | | // and set up mark as char array |
722 | | char * COIN_RESTRICT marked = reinterpret_cast<char *> (index + columnArray->capacity()); |
723 | | double * COIN_RESTRICT array2 = y->denseVector(); |
| 672 | void ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar, |
| 673 | const CoinIndexedVector *rowArray, |
| 674 | CoinIndexedVector *y, |
| 675 | CoinIndexedVector *columnArray) const |
| 676 | { |
| 677 | columnArray->clear(); |
| 678 | double *COIN_RESTRICT pi = rowArray->denseVector(); |
| 679 | int numberNonZero = 0; |
| 680 | int *COIN_RESTRICT index = columnArray->getIndices(); |
| 681 | double *COIN_RESTRICT array = columnArray->denseVector(); |
| 682 | int numberInRowArray = rowArray->getNumElements(); |
| 683 | // maybe I need one in OsiSimplex |
| 684 | double zeroTolerance = model->zeroTolerance(); |
| 685 | const int *COIN_RESTRICT column = indices_; |
| 686 | const CoinBigIndex *COIN_RESTRICT startPositive = startPositive_; |
| 687 | const CoinBigIndex *COIN_RESTRICT startNegative = startNegative_; |
| 688 | const int *COIN_RESTRICT whichRow = rowArray->getIndices(); |
| 689 | bool packed = rowArray->packedMode(); |
| 690 | if (numberInRowArray > 2) { |
| 691 | // do by rows |
| 692 | int iRow; |
| 693 | double *COIN_RESTRICT markVector = y->denseVector(); // probably empty .. but |
| 694 | int numberOriginal = 0; |
| 695 | int i; |
| 696 | if (packed) { |
| 697 | CoinBigIndex numberCovered = 0; |
| 698 | int numberColumns = getNumCols(); |
| 699 | bool sparse = true; |
| 700 | int target = 1 * numberColumns; |
| 701 | for (int i = 0; i < numberInRowArray; i++) { |
| 702 | int iRow = whichRow[i]; |
| 703 | numberCovered += startPositive[iRow + 1] - startPositive[iRow]; |
| 704 | if (numberCovered > target) { |
| 705 | sparse = false; |
| 706 | break; |
| 707 | } |
| 708 | } |
| 709 | numberNonZero = 0; |
| 710 | if (sparse) { |
| 711 | // and set up mark as char array |
| 712 | char *COIN_RESTRICT marked = reinterpret_cast< char * >(index + columnArray->capacity()); |
| 713 | double *COIN_RESTRICT array2 = y->denseVector(); |
732 | | if ((otherFlags_&1)==0||!doPlusOnes) { |
733 | | #endif |
734 | | for (int i = 0; i < numberInRowArray; i++) { |
735 | | iRow = whichRow[i]; |
736 | | double value = pi[i] * scalar; |
737 | | CoinBigIndex j; |
738 | | for (j = startPositive[iRow]; j < startNegative[iRow]; j++) { |
739 | | int iColumn = column[j]; |
740 | | if (!marked[iColumn]) { |
741 | | marked[iColumn] = 1; |
742 | | index[numberNonZero++] = iColumn; |
743 | | } |
744 | | array2[iColumn] += value; |
745 | | } |
746 | | for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) { |
747 | | int iColumn = column[j]; |
748 | | if (!marked[iColumn]) { |
749 | | marked[iColumn] = 1; |
750 | | index[numberNonZero++] = iColumn; |
751 | | } |
752 | | array2[iColumn] -= value; |
753 | | } |
754 | | } |
| 722 | if ((otherFlags_ & 1) == 0 || !doPlusOnes) { |
| 723 | #endif |
| 724 | for (int i = 0; i < numberInRowArray; i++) { |
| 725 | iRow = whichRow[i]; |
| 726 | double value = pi[i] * scalar; |
| 727 | CoinBigIndex j; |
| 728 | for (j = startPositive[iRow]; j < startNegative[iRow]; j++) { |
| 729 | int iColumn = column[j]; |
| 730 | if (!marked[iColumn]) { |
| 731 | marked[iColumn] = 1; |
| 732 | index[numberNonZero++] = iColumn; |
| 733 | } |
| 734 | array2[iColumn] += value; |
| 735 | } |
| 736 | for (j = startNegative[iRow]; j < startPositive[iRow + 1]; j++) { |
| 737 | int iColumn = column[j]; |
| 738 | if (!marked[iColumn]) { |
| 739 | marked[iColumn] = 1; |
| 740 | index[numberNonZero++] = iColumn; |
| 741 | } |
| 742 | array2[iColumn] -= value; |
| 743 | } |
| 744 | } |
756 | | } else { |
757 | | // plus one |
758 | | oneit(4); |
759 | | for (int i = 0; i < numberInRowArray; i++) { |
760 | | iRow = whichRow[i]; |
761 | | double value = pi[i] * scalar; |
762 | | CoinBigIndex j; |
763 | | for (j = startPositive[iRow]; j < startPositive[iRow+1]; j++) { |
764 | | int iColumn = column[j]; |
765 | | if (!marked[iColumn]) { |
766 | | marked[iColumn] = 1; |
767 | | index[numberNonZero++] = iColumn; |
768 | | } |
769 | | array2[iColumn] += value; |
770 | | } |
771 | | } |
772 | | } |
773 | | #endif |
774 | | // get rid of tiny values and zero out marked |
775 | | numberOriginal = numberNonZero; |
776 | | numberNonZero = 0; |
777 | | for (int i = 0; i < numberOriginal; i++) { |
778 | | int iColumn = index[i]; |
779 | | if (marked[iColumn]) { |
780 | | double value = array2[iColumn]; |
781 | | array2[iColumn] = 0.0; |
782 | | marked[iColumn] = 0; |
783 | | if (fabs(value) > zeroTolerance) { |
784 | | array[numberNonZero] = value; |
785 | | index[numberNonZero++] = iColumn; |
786 | | } |
787 | | } |
788 | | } |
789 | | } else { |
790 | | // not sparse |
| 746 | } else { |
| 747 | // plus one |
| 748 | oneit(4); |
| 749 | for (int i = 0; i < numberInRowArray; i++) { |
| 750 | iRow = whichRow[i]; |
| 751 | double value = pi[i] * scalar; |
| 752 | CoinBigIndex j; |
| 753 | for (j = startPositive[iRow]; j < startPositive[iRow + 1]; j++) { |
| 754 | int iColumn = column[j]; |
| 755 | if (!marked[iColumn]) { |
| 756 | marked[iColumn] = 1; |
| 757 | index[numberNonZero++] = iColumn; |
| 758 | } |
| 759 | array2[iColumn] += value; |
| 760 | } |
| 761 | } |
| 762 | } |
| 763 | #endif |
| 764 | // get rid of tiny values and zero out marked |
| 765 | numberOriginal = numberNonZero; |
| 766 | numberNonZero = 0; |
| 767 | for (int i = 0; i < numberOriginal; i++) { |
| 768 | int iColumn = index[i]; |
| 769 | if (marked[iColumn]) { |
| 770 | double value = array2[iColumn]; |
| 771 | array2[iColumn] = 0.0; |
| 772 | marked[iColumn] = 0; |
| 773 | if (fabs(value) > zeroTolerance) { |
| 774 | array[numberNonZero] = value; |
| 775 | index[numberNonZero++] = iColumn; |
| 776 | } |
| 777 | } |
| 778 | } |
| 779 | } else { |
| 780 | // not sparse |
808 | | } else { |
809 | | // plus one |
810 | | oneit(5); |
811 | | for (int i = 0; i < numberInRowArray; i++) { |
812 | | iRow = whichRow[i]; |
813 | | double value = pi[i] * scalar; |
814 | | CoinBigIndex j; |
815 | | for (j = startPositive[iRow]; j < startPositive[iRow+1]; j++) { |
816 | | int iColumn = column[j]; |
817 | | array[iColumn] += value; |
818 | | } |
819 | | } |
820 | | } |
821 | | #endif |
822 | | // get rid of tiny values and count |
823 | | for (int i = 0; i < numberColumns; i++) { |
824 | | double value = array[i]; |
825 | | if (value) { |
826 | | array[i] = 0.0; |
827 | | if (fabs(value) > zeroTolerance) { |
828 | | array[numberNonZero] = value; |
829 | | index[numberNonZero++] = i; |
830 | | } |
831 | | } |
832 | | } |
833 | | } |
| 798 | } else { |
| 799 | // plus one |
| 800 | oneit(5); |
| 801 | for (int i = 0; i < numberInRowArray; i++) { |
| 802 | iRow = whichRow[i]; |
| 803 | double value = pi[i] * scalar; |
| 804 | CoinBigIndex j; |
| 805 | for (j = startPositive[iRow]; j < startPositive[iRow + 1]; j++) { |
| 806 | int iColumn = column[j]; |
| 807 | array[iColumn] += value; |
| 808 | } |
| 809 | } |
| 810 | } |
| 811 | #endif |
| 812 | // get rid of tiny values and count |
| 813 | for (int i = 0; i < numberColumns; i++) { |
| 814 | double value = array[i]; |
| 815 | if (value) { |
| 816 | array[i] = 0.0; |
| 817 | if (fabs(value) > zeroTolerance) { |
| 818 | array[numberNonZero] = value; |
| 819 | index[numberNonZero++] = i; |
| 820 | } |
| 821 | } |
| 822 | } |
| 823 | } |
| 824 | } else { |
| 825 | numberNonZero = 0; |
| 826 | // and set up mark as char array |
| 827 | char *COIN_RESTRICT marked = reinterpret_cast< char * >(markVector); |
| 828 | for (i = 0; i < numberOriginal; i++) { |
| 829 | int iColumn = index[i]; |
| 830 | marked[iColumn] = 0; |
| 831 | } |
| 832 | for (i = 0; i < numberInRowArray; i++) { |
| 833 | iRow = whichRow[i]; |
| 834 | double value = pi[iRow] * scalar; |
| 835 | CoinBigIndex j; |
| 836 | for (j = startPositive[iRow]; j < startNegative[iRow]; j++) { |
| 837 | int iColumn = column[j]; |
| 838 | if (!marked[iColumn]) { |
| 839 | marked[iColumn] = 1; |
| 840 | index[numberNonZero++] = iColumn; |
| 841 | } |
| 842 | array[iColumn] += value; |
| 843 | } |
| 844 | for (j = startNegative[iRow]; j < startPositive[iRow + 1]; j++) { |
| 845 | int iColumn = column[j]; |
| 846 | if (!marked[iColumn]) { |
| 847 | marked[iColumn] = 1; |
| 848 | index[numberNonZero++] = iColumn; |
| 849 | } |
| 850 | array[iColumn] -= value; |
| 851 | } |
| 852 | } |
| 853 | // get rid of tiny values and zero out marked |
| 854 | numberOriginal = numberNonZero; |
| 855 | numberNonZero = 0; |
| 856 | for (i = 0; i < numberOriginal; i++) { |
| 857 | int iColumn = index[i]; |
| 858 | marked[iColumn] = 0; |
| 859 | if (fabs(array[iColumn]) > zeroTolerance) { |
| 860 | index[numberNonZero++] = iColumn; |
| 861 | } else { |
| 862 | array[iColumn] = 0.0; |
| 863 | } |
| 864 | } |
| 865 | } |
| 866 | } else if (numberInRowArray == 2) { |
| 867 | /* do by rows when two rows (do longer first when not packed |
| 868 | and shorter first if packed */ |
| 869 | int iRow0 = whichRow[0]; |
| 870 | int iRow1 = whichRow[1]; |
| 871 | CoinBigIndex j; |
| 872 | if (packed) { |
| 873 | double pi0 = pi[0]; |
| 874 | double pi1 = pi[1]; |
| 875 | if (startPositive[iRow0 + 1] - startPositive[iRow0] > startPositive[iRow1 + 1] - startPositive[iRow1]) { |
| 876 | int temp = iRow0; |
| 877 | iRow0 = iRow1; |
| 878 | iRow1 = temp; |
| 879 | pi0 = pi1; |
| 880 | pi1 = pi[0]; |
| 881 | } |
| 882 | // and set up mark as char array |
| 883 | char *COIN_RESTRICT marked = reinterpret_cast< char * >(index + columnArray->capacity()); |
| 884 | int *COIN_RESTRICT lookup = y->getIndices(); |
| 885 | double value = pi0 * scalar; |
| 886 | int numberOriginal; |
| 887 | #ifdef CLP_PLUS_ONE_MATRIX |
| 888 | if ((otherFlags_ & 1) == 0 || !doPlusOnes) { |
| 889 | #endif |
| 890 | for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) { |
| 891 | int iColumn = column[j]; |
| 892 | array[numberNonZero] = value; |
| 893 | marked[iColumn] = 1; |
| 894 | lookup[iColumn] = numberNonZero; |
| 895 | index[numberNonZero++] = iColumn; |
| 896 | } |
| 897 | for (j = startNegative[iRow0]; j < startPositive[iRow0 + 1]; j++) { |
| 898 | int iColumn = column[j]; |
| 899 | array[numberNonZero] = -value; |
| 900 | marked[iColumn] = 1; |
| 901 | lookup[iColumn] = numberNonZero; |
| 902 | index[numberNonZero++] = iColumn; |
| 903 | } |
| 904 | numberOriginal = numberNonZero; |
| 905 | value = pi1 * scalar; |
| 906 | for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) { |
| 907 | int iColumn = column[j]; |
| 908 | if (marked[iColumn]) { |
| 909 | int iLookup = lookup[iColumn]; |
| 910 | array[iLookup] += value; |
835 | | numberNonZero = 0; |
836 | | // and set up mark as char array |
837 | | char * COIN_RESTRICT marked = reinterpret_cast<char *> (markVector); |
838 | | for (i = 0; i < numberOriginal; i++) { |
839 | | int iColumn = index[i]; |
840 | | marked[iColumn] = 0; |
841 | | } |
842 | | for (i = 0; i < numberInRowArray; i++) { |
843 | | iRow = whichRow[i]; |
844 | | double value = pi[iRow] * scalar; |
845 | | CoinBigIndex j; |
846 | | for (j = startPositive[iRow]; j < startNegative[iRow]; j++) { |
847 | | int iColumn = column[j]; |
848 | | if (!marked[iColumn]) { |
849 | | marked[iColumn] = 1; |
850 | | index[numberNonZero++] = iColumn; |
851 | | } |
852 | | array[iColumn] += value; |
853 | | } |
854 | | for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) { |
855 | | int iColumn = column[j]; |
856 | | if (!marked[iColumn]) { |
857 | | marked[iColumn] = 1; |
858 | | index[numberNonZero++] = iColumn; |
859 | | } |
860 | | array[iColumn] -= value; |
861 | | } |
862 | | } |
863 | | // get rid of tiny values and zero out marked |
864 | | numberOriginal = numberNonZero; |
865 | | numberNonZero = 0; |
866 | | for (i = 0; i < numberOriginal; i++) { |
867 | | int iColumn = index[i]; |
868 | | marked[iColumn] = 0; |
869 | | if (fabs(array[iColumn]) > zeroTolerance) { |
870 | | index[numberNonZero++] = iColumn; |
871 | | } else { |
872 | | array[iColumn] = 0.0; |
873 | | } |
874 | | } |
| 912 | if (fabs(value) > zeroTolerance) { |
| 913 | array[numberNonZero] = value; |
| 914 | index[numberNonZero++] = iColumn; |
| 915 | } |
899 | | if ((otherFlags_&1)==0||!doPlusOnes) { |
900 | | #endif |
901 | | for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) { |
902 | | int iColumn = column[j]; |
903 | | array[numberNonZero] = value; |
904 | | marked[iColumn] = 1; |
905 | | lookup[iColumn] = numberNonZero; |
906 | | index[numberNonZero++] = iColumn; |
907 | | } |
908 | | for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) { |
909 | | int iColumn = column[j]; |
910 | | array[numberNonZero] = -value; |
911 | | marked[iColumn] = 1; |
912 | | lookup[iColumn] = numberNonZero; |
913 | | index[numberNonZero++] = iColumn; |
914 | | } |
915 | | numberOriginal = numberNonZero; |
916 | | value = pi1 * scalar; |
917 | | for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) { |
918 | | int iColumn = column[j]; |
919 | | if (marked[iColumn]) { |
920 | | int iLookup = lookup[iColumn]; |
921 | | array[iLookup] += value; |
922 | | } else { |
923 | | if (fabs(value) > zeroTolerance) { |
924 | | array[numberNonZero] = value; |
925 | | index[numberNonZero++] = iColumn; |
926 | | } |
927 | | } |
928 | | } |
929 | | for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) { |
930 | | int iColumn = column[j]; |
931 | | if (marked[iColumn]) { |
932 | | int iLookup = lookup[iColumn]; |
933 | | array[iLookup] -= value; |
934 | | } else { |
935 | | if (fabs(value) > zeroTolerance) { |
936 | | array[numberNonZero] = -value; |
937 | | index[numberNonZero++] = iColumn; |
938 | | } |
939 | | } |
940 | | } |
| 931 | } else { |
| 932 | // plus one |
| 933 | oneit(7); |
| 934 | for (j = startPositive[iRow0]; j < startPositive[iRow0 + 1]; j++) { |
| 935 | int iColumn = column[j]; |
| 936 | array[numberNonZero] = value; |
| 937 | marked[iColumn] = 1; |
| 938 | lookup[iColumn] = numberNonZero; |
| 939 | index[numberNonZero++] = iColumn; |
| 940 | } |
| 941 | numberOriginal = numberNonZero; |
| 942 | value = pi1 * scalar; |
| 943 | for (j = startPositive[iRow1]; j < startPositive[iRow1 + 1]; j++) { |
| 944 | int iColumn = column[j]; |
| 945 | if (marked[iColumn]) { |
| 946 | int iLookup = lookup[iColumn]; |
| 947 | array[iLookup] += value; |
| 948 | } else { |
| 949 | if (fabs(value) > zeroTolerance) { |
| 950 | array[numberNonZero] = value; |
| 951 | index[numberNonZero++] = iColumn; |
| 952 | } |
| 953 | } |
| 954 | } |
| 955 | } |
| 956 | #endif |
| 957 | // get rid of tiny values and zero out marked |
| 958 | int nDelete = 0; |
| 959 | for (j = 0; j < numberOriginal; j++) { |
| 960 | int iColumn = index[j]; |
| 961 | marked[iColumn] = 0; |
| 962 | if (fabs(array[j]) <= zeroTolerance) |
| 963 | nDelete++; |
| 964 | } |
| 965 | if (nDelete) { |
| 966 | numberOriginal = numberNonZero; |
| 967 | numberNonZero = 0; |
| 968 | for (j = 0; j < numberOriginal; j++) { |
| 969 | int iColumn = index[j]; |
| 970 | double value = array[j]; |
| 971 | array[j] = 0.0; |
| 972 | if (fabs(value) > zeroTolerance) { |
| 973 | array[numberNonZero] = value; |
| 974 | index[numberNonZero++] = iColumn; |
| 975 | } |
| 976 | } |
| 977 | } |
| 978 | } else { |
| 979 | if (startPositive[iRow0 + 1] - startPositive[iRow0] < startPositive[iRow1 + 1] - startPositive[iRow1]) { |
| 980 | int temp = iRow0; |
| 981 | iRow0 = iRow1; |
| 982 | iRow1 = temp; |
| 983 | } |
| 984 | int numberOriginal; |
| 985 | int i; |
| 986 | numberNonZero = 0; |
| 987 | double value; |
| 988 | value = pi[iRow0] * scalar; |
| 989 | CoinBigIndex j; |
| 990 | for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) { |
| 991 | int iColumn = column[j]; |
| 992 | index[numberNonZero++] = iColumn; |
| 993 | array[iColumn] = value; |
| 994 | } |
| 995 | for (j = startNegative[iRow0]; j < startPositive[iRow0 + 1]; j++) { |
| 996 | int iColumn = column[j]; |
| 997 | index[numberNonZero++] = iColumn; |
| 998 | array[iColumn] = -value; |
| 999 | } |
| 1000 | value = pi[iRow1] * scalar; |
| 1001 | for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) { |
| 1002 | int iColumn = column[j]; |
| 1003 | double value2 = array[iColumn]; |
| 1004 | if (value2) { |
| 1005 | value2 += value; |
| 1006 | } else { |
| 1007 | value2 = value; |
| 1008 | index[numberNonZero++] = iColumn; |
| 1009 | } |
| 1010 | array[iColumn] = value2; |
| 1011 | } |
| 1012 | for (j = startNegative[iRow1]; j < startPositive[iRow1 + 1]; j++) { |
| 1013 | int iColumn = column[j]; |
| 1014 | double value2 = array[iColumn]; |
| 1015 | if (value2) { |
| 1016 | value2 -= value; |
| 1017 | } else { |
| 1018 | value2 = -value; |
| 1019 | index[numberNonZero++] = iColumn; |
| 1020 | } |
| 1021 | array[iColumn] = value2; |
| 1022 | } |
| 1023 | // get rid of tiny values and zero out marked |
| 1024 | numberOriginal = numberNonZero; |
| 1025 | numberNonZero = 0; |
| 1026 | for (i = 0; i < numberOriginal; i++) { |
| 1027 | int iColumn = index[i]; |
| 1028 | if (fabs(array[iColumn]) > zeroTolerance) { |
| 1029 | index[numberNonZero++] = iColumn; |
| 1030 | } else { |
| 1031 | array[iColumn] = 0.0; |
| 1032 | } |
| 1033 | } |
| 1034 | } |
| 1035 | } else if (numberInRowArray == 1) { |
| 1036 | // Just one row |
| 1037 | int iRow = rowArray->getIndices()[0]; |
| 1038 | numberNonZero = 0; |
| 1039 | double value; |
| 1040 | iRow = whichRow[0]; |
| 1041 | CoinBigIndex j; |
| 1042 | if (packed) { |
| 1043 | value = pi[0] * scalar; |
| 1044 | if (fabs(value) > zeroTolerance) { |
942 | | } else { |
943 | | // plus one |
944 | | oneit(7); |
945 | | for (j = startPositive[iRow0]; j < startPositive[iRow0+1]; j++) { |
946 | | int iColumn = column[j]; |
947 | | array[numberNonZero] = value; |
948 | | marked[iColumn] = 1; |
949 | | lookup[iColumn] = numberNonZero; |
950 | | index[numberNonZero++] = iColumn; |
951 | | } |
952 | | numberOriginal = numberNonZero; |
953 | | value = pi1 * scalar; |
954 | | for (j = startPositive[iRow1]; j < startPositive[iRow1+1]; j++) { |
955 | | int iColumn = column[j]; |
956 | | if (marked[iColumn]) { |
957 | | int iLookup = lookup[iColumn]; |
958 | | array[iLookup] += value; |
959 | | } else { |
960 | | if (fabs(value) > zeroTolerance) { |
961 | | array[numberNonZero] = value; |
962 | | index[numberNonZero++] = iColumn; |
963 | | } |
964 | | } |
965 | | } |
966 | | } |
967 | | #endif |
968 | | // get rid of tiny values and zero out marked |
969 | | int nDelete = 0; |
970 | | for (j = 0; j < numberOriginal; j++) { |
971 | | int iColumn = index[j]; |
972 | | marked[iColumn] = 0; |
973 | | if (fabs(array[j]) <= zeroTolerance) |
974 | | nDelete++; |
975 | | } |
976 | | if (nDelete) { |
977 | | numberOriginal = numberNonZero; |
978 | | numberNonZero = 0; |
979 | | for (j = 0; j < numberOriginal; j++) { |
980 | | int iColumn = index[j]; |
981 | | double value = array[j]; |
982 | | array[j] = 0.0; |
983 | | if (fabs(value) > zeroTolerance) { |
984 | | array[numberNonZero] = value; |
985 | | index[numberNonZero++] = iColumn; |
986 | | } |
987 | | } |
988 | | } |
989 | | } else { |
990 | | if (startPositive[iRow0+1] - startPositive[iRow0] < |
991 | | startPositive[iRow1+1] - startPositive[iRow1]) { |
992 | | int temp = iRow0; |
993 | | iRow0 = iRow1; |
994 | | iRow1 = temp; |
995 | | } |
996 | | int numberOriginal; |
997 | | int i; |
998 | | numberNonZero = 0; |
999 | | double value; |
1000 | | value = pi[iRow0] * scalar; |
1001 | | CoinBigIndex j; |
1002 | | for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) { |
1003 | | int iColumn = column[j]; |
1004 | | index[numberNonZero++] = iColumn; |
1005 | | array[iColumn] = value; |
1006 | | } |
1007 | | for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) { |
1008 | | int iColumn = column[j]; |
1009 | | index[numberNonZero++] = iColumn; |
1010 | | array[iColumn] = -value; |
1011 | | } |
1012 | | value = pi[iRow1] * scalar; |
1013 | | for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) { |
1014 | | int iColumn = column[j]; |
1015 | | double value2 = array[iColumn]; |
1016 | | if (value2) { |
1017 | | value2 += value; |
1018 | | } else { |
1019 | | value2 = value; |
1020 | | index[numberNonZero++] = iColumn; |
1021 | | } |
1022 | | array[iColumn] = value2; |
1023 | | } |
1024 | | for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) { |
1025 | | int iColumn = column[j]; |
1026 | | double value2 = array[iColumn]; |
1027 | | if (value2) { |
1028 | | value2 -= value; |
1029 | | } else { |
1030 | | value2 = -value; |
1031 | | index[numberNonZero++] = iColumn; |
1032 | | } |
1033 | | array[iColumn] = value2; |
1034 | | } |
1035 | | // get rid of tiny values and zero out marked |
1036 | | numberOriginal = numberNonZero; |
1037 | | numberNonZero = 0; |
1038 | | for (i = 0; i < numberOriginal; i++) { |
1039 | | int iColumn = index[i]; |
1040 | | if (fabs(array[iColumn]) > zeroTolerance) { |
1041 | | index[numberNonZero++] = iColumn; |
1042 | | } else { |
1043 | | array[iColumn] = 0.0; |
1044 | | } |
1045 | | } |
| 1046 | if ((otherFlags_ & 1) == 0 || !doPlusOnes) { |
| 1047 | #endif |
| 1048 | for (j = startPositive[iRow]; j < startNegative[iRow]; j++) { |
| 1049 | int iColumn = column[j]; |
| 1050 | array[numberNonZero] = value; |
| 1051 | index[numberNonZero++] = iColumn; |
1058 | | if ((otherFlags_&1)==0||!doPlusOnes) { |
1059 | | #endif |
1060 | | for (j = startPositive[iRow]; j < startNegative[iRow]; j++) { |
1061 | | int iColumn = column[j]; |
1062 | | array[numberNonZero] = value; |
1063 | | index[numberNonZero++] = iColumn; |
1064 | | } |
1065 | | for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) { |
1066 | | int iColumn = column[j]; |
1067 | | array[numberNonZero] = -value; |
1068 | | index[numberNonZero++] = iColumn; |
1069 | | } |
1070 | | #ifdef CLP_PLUS_ONE_MATRIX |
1071 | | } else { |
1072 | | // plus one |
1073 | | oneit(9); |
1074 | | for (j = startPositive[iRow]; j < startPositive[iRow+1]; j++) { |
1075 | | int iColumn = column[j]; |
1076 | | array[numberNonZero] = value; |
1077 | | index[numberNonZero++] = iColumn; |
1078 | | } |
1079 | | } |
1080 | | #endif |
1081 | | } |
1082 | | } else { |
1083 | | value = pi[iRow] * scalar; |
1084 | | if (fabs(value) > zeroTolerance) { |
1085 | | for (j = startPositive[iRow]; j < startNegative[iRow]; j++) { |
1086 | | int iColumn = column[j]; |
1087 | | array[iColumn] = value; |
1088 | | index[numberNonZero++] = iColumn; |
1089 | | } |
1090 | | for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) { |
1091 | | int iColumn = column[j]; |
1092 | | array[iColumn] = -value; |
1093 | | index[numberNonZero++] = iColumn; |
1094 | | } |
1095 | | } |
| 1059 | } else { |
| 1060 | // plus one |
| 1061 | oneit(9); |
| 1062 | for (j = startPositive[iRow]; j < startPositive[iRow + 1]; j++) { |
| 1063 | int iColumn = column[j]; |
| 1064 | array[numberNonZero] = value; |
| 1065 | index[numberNonZero++] = iColumn; |
1105 | | void |
1106 | | ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex * , |
1107 | | const CoinIndexedVector * rowArray, |
1108 | | const CoinIndexedVector * y, |
1109 | | CoinIndexedVector * columnArray) const |
1110 | | { |
1111 | | columnArray->clear(); |
1112 | | double * COIN_RESTRICT pi = rowArray->denseVector(); |
1113 | | double * COIN_RESTRICT array = columnArray->denseVector(); |
1114 | | int jColumn; |
1115 | | int numberToDo = y->getNumElements(); |
1116 | | const int * COIN_RESTRICT which = y->getIndices(); |
1117 | | assert (!rowArray->packedMode()); |
1118 | | columnArray->setPacked(); |
| 1093 | void ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex *, |
| 1094 | const CoinIndexedVector *rowArray, |
| 1095 | const CoinIndexedVector *y, |
| 1096 | CoinIndexedVector *columnArray) const |
| 1097 | { |
| 1098 | columnArray->clear(); |
| 1099 | double *COIN_RESTRICT pi = rowArray->denseVector(); |
| 1100 | double *COIN_RESTRICT array = columnArray->denseVector(); |
| 1101 | int jColumn; |
| 1102 | int numberToDo = y->getNumElements(); |
| 1103 | const int *COIN_RESTRICT which = y->getIndices(); |
| 1104 | assert(!rowArray->packedMode()); |
| 1105 | columnArray->setPacked(); |
1165 | | if (numberElements>COIN_INT_MAX) { |
1166 | | printf("Factorization too large\n"); |
1167 | | abort(); |
1168 | | } |
1169 | | #endif |
1170 | | return static_cast<int>(numberElements); |
1171 | | } |
1172 | | void |
1173 | | ClpPlusMinusOneMatrix::fillBasis(ClpSimplex * , |
1174 | | const int * whichColumn, |
1175 | | int & numberColumnBasic, |
1176 | | int * indexRowU, int * start, |
1177 | | int * rowCount, int * columnCount, |
1178 | | CoinFactorizationDouble * elementU) |
1179 | | { |
1180 | | int i; |
1181 | | CoinBigIndex numberElements = start[0]; |
1182 | | assert (columnOrdered_); |
1183 | | for (i = 0; i < numberColumnBasic; i++) { |
1184 | | int iColumn = whichColumn[i]; |
1185 | | CoinBigIndex j = startPositive_[iColumn]; |
1186 | | for (; j < startNegative_[iColumn]; j++) { |
1187 | | int iRow = indices_[j]; |
1188 | | indexRowU[numberElements] = iRow; |
1189 | | rowCount[iRow]++; |
1190 | | elementU[numberElements++] = 1.0; |
1191 | | } |
1192 | | for (; j < startPositive_[iColumn+1]; j++) { |
1193 | | int iRow = indices_[j]; |
1194 | | indexRowU[numberElements] = iRow; |
1195 | | rowCount[iRow]++; |
1196 | | elementU[numberElements++] = -1.0; |
1197 | | } |
1198 | | start[i+1] = static_cast<int>(numberElements); |
1199 | | columnCount[i] = static_cast<int>(numberElements - start[i]); |
1200 | | } |
| 1151 | if (numberElements > COIN_INT_MAX) { |
| 1152 | printf("Factorization too large\n"); |
| 1153 | abort(); |
| 1154 | } |
| 1155 | #endif |
| 1156 | return static_cast< int >(numberElements); |
| 1157 | } |
| 1158 | void ClpPlusMinusOneMatrix::fillBasis(ClpSimplex *, |
| 1159 | const int *whichColumn, |
| 1160 | int &numberColumnBasic, |
| 1161 | int *indexRowU, int *start, |
| 1162 | int *rowCount, int *columnCount, |
| 1163 | CoinFactorizationDouble *elementU) |
| 1164 | { |
| 1165 | int i; |
| 1166 | CoinBigIndex numberElements = start[0]; |
| 1167 | assert(columnOrdered_); |
| 1168 | for (i = 0; i < numberColumnBasic; i++) { |
| 1169 | int iColumn = whichColumn[i]; |
| 1170 | CoinBigIndex j = startPositive_[iColumn]; |
| 1171 | for (; j < startNegative_[iColumn]; j++) { |
| 1172 | int iRow = indices_[j]; |
| 1173 | indexRowU[numberElements] = iRow; |
| 1174 | rowCount[iRow]++; |
| 1175 | elementU[numberElements++] = 1.0; |
| 1176 | } |
| 1177 | for (; j < startPositive_[iColumn + 1]; j++) { |
| 1178 | int iRow = indices_[j]; |
| 1179 | indexRowU[numberElements] = iRow; |
| 1180 | rowCount[iRow]++; |
| 1181 | elementU[numberElements++] = -1.0; |
| 1182 | } |
| 1183 | start[i + 1] = static_cast< int >(numberElements); |
| 1184 | columnCount[i] = static_cast< int >(numberElements - start[i]); |
| 1185 | } |
1208 | | void |
1209 | | ClpPlusMinusOneMatrix::unpack(const ClpSimplex * , |
1210 | | CoinIndexedVector * rowArray, |
1211 | | int iColumn) const |
1212 | | { |
1213 | | CoinBigIndex j = startPositive_[iColumn]; |
1214 | | for (; j < startNegative_[iColumn]; j++) { |
1215 | | int iRow = indices_[j]; |
1216 | | rowArray->add(iRow, 1.0); |
1217 | | } |
1218 | | for (; j < startPositive_[iColumn+1]; j++) { |
1219 | | int iRow = indices_[j]; |
1220 | | rowArray->add(iRow, -1.0); |
1221 | | } |
| 1193 | void ClpPlusMinusOneMatrix::unpack(const ClpSimplex *, |
| 1194 | CoinIndexedVector *rowArray, |
| 1195 | int iColumn) const |
| 1196 | { |
| 1197 | CoinBigIndex j = startPositive_[iColumn]; |
| 1198 | for (; j < startNegative_[iColumn]; j++) { |
| 1199 | int iRow = indices_[j]; |
| 1200 | rowArray->add(iRow, 1.0); |
| 1201 | } |
| 1202 | for (; j < startPositive_[iColumn + 1]; j++) { |
| 1203 | int iRow = indices_[j]; |
| 1204 | rowArray->add(iRow, -1.0); |
| 1205 | } |
1227 | | void |
1228 | | ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * , |
1229 | | CoinIndexedVector * rowArray, |
1230 | | int iColumn) const |
1231 | | { |
1232 | | int * COIN_RESTRICT index = rowArray->getIndices(); |
1233 | | double * COIN_RESTRICT array = rowArray->denseVector(); |
1234 | | int number = 0; |
1235 | | CoinBigIndex j = startPositive_[iColumn]; |
1236 | | for (; j < startNegative_[iColumn]; j++) { |
1237 | | int iRow = indices_[j]; |
1238 | | array[number] = 1.0; |
1239 | | index[number++] = iRow; |
1240 | | } |
1241 | | for (; j < startPositive_[iColumn+1]; j++) { |
1242 | | int iRow = indices_[j]; |
1243 | | array[number] = -1.0; |
1244 | | index[number++] = iRow; |
1245 | | } |
1246 | | rowArray->setNumElements(number); |
1247 | | rowArray->setPackedMode(true); |
| 1211 | void ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex *, |
| 1212 | CoinIndexedVector *rowArray, |
| 1213 | int iColumn) const |
| 1214 | { |
| 1215 | int *COIN_RESTRICT index = rowArray->getIndices(); |
| 1216 | double *COIN_RESTRICT array = rowArray->denseVector(); |
| 1217 | int number = 0; |
| 1218 | CoinBigIndex j = startPositive_[iColumn]; |
| 1219 | for (; j < startNegative_[iColumn]; j++) { |
| 1220 | int iRow = indices_[j]; |
| 1221 | array[number] = 1.0; |
| 1222 | index[number++] = iRow; |
| 1223 | } |
| 1224 | for (; j < startPositive_[iColumn + 1]; j++) { |
| 1225 | int iRow = indices_[j]; |
| 1226 | array[number] = -1.0; |
| 1227 | index[number++] = iRow; |
| 1228 | } |
| 1229 | rowArray->setNumElements(number); |
| 1230 | rowArray->setPackedMode(true); |
1251 | | void |
1252 | | ClpPlusMinusOneMatrix::add(const ClpSimplex * , CoinIndexedVector * rowArray, |
1253 | | int iColumn, double multiplier) const |
1254 | | { |
1255 | | CoinBigIndex j = startPositive_[iColumn]; |
1256 | | for (; j < startNegative_[iColumn]; j++) { |
1257 | | int iRow = indices_[j]; |
1258 | | rowArray->quickAdd(iRow, multiplier); |
1259 | | } |
1260 | | for (; j < startPositive_[iColumn+1]; j++) { |
1261 | | int iRow = indices_[j]; |
1262 | | rowArray->quickAdd(iRow, -multiplier); |
1263 | | } |
| 1234 | void ClpPlusMinusOneMatrix::add(const ClpSimplex *, CoinIndexedVector *rowArray, |
| 1235 | int iColumn, double multiplier) const |
| 1236 | { |
| 1237 | CoinBigIndex j = startPositive_[iColumn]; |
| 1238 | for (; j < startNegative_[iColumn]; j++) { |
| 1239 | int iRow = indices_[j]; |
| 1240 | rowArray->quickAdd(iRow, multiplier); |
| 1241 | } |
| 1242 | for (; j < startPositive_[iColumn + 1]; j++) { |
| 1243 | int iRow = indices_[j]; |
| 1244 | rowArray->quickAdd(iRow, -multiplier); |
| 1245 | } |
1285 | | if (!matrix_) { |
1286 | | int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_; |
1287 | | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
1288 | | CoinBigIndex numberElements = startPositive_[numberMajor]; |
1289 | | double * elements = new double [numberElements]; |
1290 | | CoinBigIndex j = 0; |
1291 | | int i; |
1292 | | for (i = 0; i < numberMajor; i++) { |
1293 | | for (; j < startNegative_[i]; j++) { |
1294 | | elements[j] = 1.0; |
1295 | | } |
1296 | | for (; j < startPositive_[i+1]; j++) { |
1297 | | elements[j] = -1.0; |
1298 | | } |
1299 | | } |
1300 | | matrix_ = new CoinPackedMatrix(columnOrdered_, numberMinor, numberMajor, |
1301 | | getNumElements(), |
1302 | | elements, indices_, |
1303 | | startPositive_, getVectorLengths()); |
1304 | | delete [] elements; |
1305 | | delete [] lengths_; |
1306 | | lengths_ = NULL; |
1307 | | } |
1308 | | return matrix_; |
| 1266 | if (!matrix_) { |
| 1267 | int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_; |
| 1268 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
| 1269 | CoinBigIndex numberElements = startPositive_[numberMajor]; |
| 1270 | double *elements = new double[numberElements]; |
| 1271 | CoinBigIndex j = 0; |
| 1272 | int i; |
| 1273 | for (i = 0; i < numberMajor; i++) { |
| 1274 | for (; j < startNegative_[i]; j++) { |
| 1275 | elements[j] = 1.0; |
| 1276 | } |
| 1277 | for (; j < startPositive_[i + 1]; j++) { |
| 1278 | elements[j] = -1.0; |
| 1279 | } |
| 1280 | } |
| 1281 | matrix_ = new CoinPackedMatrix(columnOrdered_, numberMinor, numberMajor, |
| 1282 | getNumElements(), |
| 1283 | elements, indices_, |
| 1284 | startPositive_, getVectorLengths()); |
| 1285 | delete[] elements; |
| 1286 | delete[] lengths_; |
| 1287 | lengths_ = NULL; |
| 1288 | } |
| 1289 | return matrix_; |
1342 | | void |
1343 | | ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel) |
1344 | | { |
1345 | | int iColumn; |
1346 | | CoinBigIndex newSize = startPositive_[numberColumns_];; |
1347 | | int numberBad = 0; |
1348 | | // Use array to make sure we can have duplicates |
1349 | | int * which = new int[numberColumns_]; |
1350 | | memset(which, 0, numberColumns_ * sizeof(int)); |
1351 | | int nDuplicate = 0; |
1352 | | for (iColumn = 0; iColumn < numDel; iColumn++) { |
1353 | | int jColumn = indDel[iColumn]; |
1354 | | if (jColumn < 0 || jColumn >= numberColumns_) { |
1355 | | numberBad++; |
1356 | | } else { |
1357 | | newSize -= startPositive_[jColumn+1] - startPositive_[jColumn]; |
1358 | | if (which[jColumn]) |
1359 | | nDuplicate++; |
1360 | | else |
1361 | | which[jColumn] = 1; |
1362 | | } |
1363 | | } |
1364 | | if (numberBad) |
1365 | | throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix"); |
1366 | | int newNumber = numberColumns_ - numDel + nDuplicate; |
1367 | | // Get rid of temporary arrays |
1368 | | delete [] lengths_; |
1369 | | lengths_ = NULL; |
1370 | | delete matrix_; |
1371 | | matrix_ = NULL; |
1372 | | CoinBigIndex * newPositive = new CoinBigIndex [newNumber+1]; |
1373 | | CoinBigIndex * newNegative = new CoinBigIndex [newNumber]; |
1374 | | int * newIndices = new int [newSize]; |
1375 | | newNumber = 0; |
1376 | | newSize = 0; |
1377 | | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
1378 | | if (!which[iColumn]) { |
1379 | | CoinBigIndex start, end; |
1380 | | CoinBigIndex i; |
1381 | | start = startPositive_[iColumn]; |
1382 | | end = startNegative_[iColumn]; |
1383 | | newPositive[newNumber] = newSize; |
1384 | | for (i = start; i < end; i++) |
1385 | | newIndices[newSize++] = indices_[i]; |
1386 | | start = startNegative_[iColumn]; |
1387 | | end = startPositive_[iColumn+1]; |
1388 | | newNegative[newNumber++] = newSize; |
1389 | | for (i = start; i < end; i++) |
1390 | | newIndices[newSize++] = indices_[i]; |
1391 | | } |
1392 | | } |
1393 | | newPositive[newNumber] = newSize; |
1394 | | delete [] which; |
1395 | | delete [] startPositive_; |
1396 | | startPositive_ = newPositive; |
1397 | | delete [] startNegative_; |
1398 | | startNegative_ = newNegative; |
1399 | | delete [] indices_; |
1400 | | indices_ = newIndices; |
1401 | | numberColumns_ = newNumber; |
| 1323 | void ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int *indDel) |
| 1324 | { |
| 1325 | int iColumn; |
| 1326 | CoinBigIndex newSize = startPositive_[numberColumns_]; |
| 1327 | ; |
| 1328 | int numberBad = 0; |
| 1329 | // Use array to make sure we can have duplicates |
| 1330 | int *which = new int[numberColumns_]; |
| 1331 | memset(which, 0, numberColumns_ * sizeof(int)); |
| 1332 | int nDuplicate = 0; |
| 1333 | for (iColumn = 0; iColumn < numDel; iColumn++) { |
| 1334 | int jColumn = indDel[iColumn]; |
| 1335 | if (jColumn < 0 || jColumn >= numberColumns_) { |
| 1336 | numberBad++; |
| 1337 | } else { |
| 1338 | newSize -= startPositive_[jColumn + 1] - startPositive_[jColumn]; |
| 1339 | if (which[jColumn]) |
| 1340 | nDuplicate++; |
| 1341 | else |
| 1342 | which[jColumn] = 1; |
| 1343 | } |
| 1344 | } |
| 1345 | if (numberBad) |
| 1346 | throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix"); |
| 1347 | int newNumber = numberColumns_ - numDel + nDuplicate; |
| 1348 | // Get rid of temporary arrays |
| 1349 | delete[] lengths_; |
| 1350 | lengths_ = NULL; |
| 1351 | delete matrix_; |
| 1352 | matrix_ = NULL; |
| 1353 | CoinBigIndex *newPositive = new CoinBigIndex[newNumber + 1]; |
| 1354 | CoinBigIndex *newNegative = new CoinBigIndex[newNumber]; |
| 1355 | int *newIndices = new int[newSize]; |
| 1356 | newNumber = 0; |
| 1357 | newSize = 0; |
| 1358 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
| 1359 | if (!which[iColumn]) { |
| 1360 | CoinBigIndex start, end; |
| 1361 | CoinBigIndex i; |
| 1362 | start = startPositive_[iColumn]; |
| 1363 | end = startNegative_[iColumn]; |
| 1364 | newPositive[newNumber] = newSize; |
| 1365 | for (i = start; i < end; i++) |
| 1366 | newIndices[newSize++] = indices_[i]; |
| 1367 | start = startNegative_[iColumn]; |
| 1368 | end = startPositive_[iColumn + 1]; |
| 1369 | newNegative[newNumber++] = newSize; |
| 1370 | for (i = start; i < end; i++) |
| 1371 | newIndices[newSize++] = indices_[i]; |
| 1372 | } |
| 1373 | } |
| 1374 | newPositive[newNumber] = newSize; |
| 1375 | delete[] which; |
| 1376 | delete[] startPositive_; |
| 1377 | startPositive_ = newPositive; |
| 1378 | delete[] startNegative_; |
| 1379 | startNegative_ = newNegative; |
| 1380 | delete[] indices_; |
| 1381 | indices_ = newIndices; |
| 1382 | numberColumns_ = newNumber; |
1404 | | void |
1405 | | ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel) |
1406 | | { |
1407 | | int iRow; |
1408 | | int numberBad = 0; |
1409 | | // Use array to make sure we can have duplicates |
1410 | | int * which = new int[numberRows_]; |
1411 | | memset(which, 0, numberRows_ * sizeof(int)); |
1412 | | int nDuplicate = 0; |
1413 | | for (iRow = 0; iRow < numDel; iRow++) { |
1414 | | int jRow = indDel[iRow]; |
1415 | | if (jRow < 0 || jRow >= numberRows_) { |
1416 | | numberBad++; |
1417 | | } else { |
1418 | | if (which[jRow]) |
1419 | | nDuplicate++; |
1420 | | else |
1421 | | which[jRow] = 1; |
1422 | | } |
1423 | | } |
1424 | | if (numberBad) |
1425 | | throw CoinError("Indices out of range", "deleteRows", "ClpPlusMinusOneMatrix"); |
1426 | | CoinBigIndex iElement; |
1427 | | CoinBigIndex numberElements = startPositive_[numberColumns_]; |
1428 | | CoinBigIndex newSize = 0; |
1429 | | for (iElement = 0; iElement < numberElements; iElement++) { |
1430 | | iRow = indices_[iElement]; |
1431 | | if (!which[iRow]) |
1432 | | newSize++; |
1433 | | } |
1434 | | int newNumber = numberRows_ - numDel + nDuplicate; |
1435 | | // Get rid of temporary arrays |
1436 | | delete [] lengths_; |
1437 | | lengths_ = NULL; |
1438 | | delete matrix_; |
1439 | | matrix_ = NULL; |
1440 | | // redo which |
1441 | | int numberRows=0; |
1442 | | for (iRow = 0; iRow < numberRows_; iRow++) { |
1443 | | if (which[iRow]) { |
1444 | | which[iRow]=-1; // not wanted |
1445 | | } else { |
1446 | | which[iRow]=numberRows; |
1447 | | numberRows++; |
1448 | | } |
1449 | | } |
1450 | | int * newIndices = new int [newSize]; |
1451 | | newSize = 0; |
1452 | | int iColumn; |
1453 | | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
1454 | | CoinBigIndex start, end; |
1455 | | CoinBigIndex i; |
1456 | | start = startPositive_[iColumn]; |
1457 | | end = startNegative_[iColumn]; |
1458 | | startPositive_[newNumber] = newSize; |
1459 | | for (i = start; i < end; i++) { |
1460 | | iRow = which[indices_[i]]; |
1461 | | if (iRow>=0) |
1462 | | newIndices[newSize++] = iRow; |
1463 | | } |
1464 | | start = startNegative_[iColumn]; |
1465 | | end = startPositive_[iColumn+1]; |
1466 | | startNegative_[newNumber] = newSize; |
1467 | | for (i = start; i < end; i++) { |
1468 | | iRow = which[indices_[i]]; |
1469 | | if (iRow>=0) |
1470 | | newIndices[newSize++] = iRow; |
1471 | | } |
1472 | | } |
1473 | | startPositive_[numberColumns_] = newSize; |
1474 | | delete [] which; |
1475 | | delete [] indices_; |
1476 | | indices_ = newIndices; |
1477 | | numberRows_ = newNumber; |
1478 | | } |
1479 | | bool |
1480 | | ClpPlusMinusOneMatrix::isColOrdered() const |
1481 | | { |
1482 | | return columnOrdered_; |
| 1385 | void ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int *indDel) |
| 1386 | { |
| 1387 | int iRow; |
| 1388 | int numberBad = 0; |
| 1389 | // Use array to make sure we can have duplicates |
| 1390 | int *which = new int[numberRows_]; |
| 1391 | memset(which, 0, numberRows_ * sizeof(int)); |
| 1392 | int nDuplicate = 0; |
| 1393 | for (iRow = 0; iRow < numDel; iRow++) { |
| 1394 | int jRow = indDel[iRow]; |
| 1395 | if (jRow < 0 || jRow >= numberRows_) { |
| 1396 | numberBad++; |
| 1397 | } else { |
| 1398 | if (which[jRow]) |
| 1399 | nDuplicate++; |
| 1400 | else |
| 1401 | which[jRow] = 1; |
| 1402 | } |
| 1403 | } |
| 1404 | if (numberBad) |
| 1405 | throw CoinError("Indices out of range", "deleteRows", "ClpPlusMinusOneMatrix"); |
| 1406 | CoinBigIndex iElement; |
| 1407 | CoinBigIndex numberElements = startPositive_[numberColumns_]; |
| 1408 | CoinBigIndex newSize = 0; |
| 1409 | for (iElement = 0; iElement < numberElements; iElement++) { |
| 1410 | iRow = indices_[iElement]; |
| 1411 | if (!which[iRow]) |
| 1412 | newSize++; |
| 1413 | } |
| 1414 | int newNumber = numberRows_ - numDel + nDuplicate; |
| 1415 | // Get rid of temporary arrays |
| 1416 | delete[] lengths_; |
| 1417 | lengths_ = NULL; |
| 1418 | delete matrix_; |
| 1419 | matrix_ = NULL; |
| 1420 | // redo which |
| 1421 | int numberRows = 0; |
| 1422 | for (iRow = 0; iRow < numberRows_; iRow++) { |
| 1423 | if (which[iRow]) { |
| 1424 | which[iRow] = -1; // not wanted |
| 1425 | } else { |
| 1426 | which[iRow] = numberRows; |
| 1427 | numberRows++; |
| 1428 | } |
| 1429 | } |
| 1430 | int *newIndices = new int[newSize]; |
| 1431 | newSize = 0; |
| 1432 | int iColumn; |
| 1433 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
| 1434 | CoinBigIndex start, end; |
| 1435 | CoinBigIndex i; |
| 1436 | start = startPositive_[iColumn]; |
| 1437 | end = startNegative_[iColumn]; |
| 1438 | startPositive_[newNumber] = newSize; |
| 1439 | for (i = start; i < end; i++) { |
| 1440 | iRow = which[indices_[i]]; |
| 1441 | if (iRow >= 0) |
| 1442 | newIndices[newSize++] = iRow; |
| 1443 | } |
| 1444 | start = startNegative_[iColumn]; |
| 1445 | end = startPositive_[iColumn + 1]; |
| 1446 | startNegative_[newNumber] = newSize; |
| 1447 | for (i = start; i < end; i++) { |
| 1448 | iRow = which[indices_[i]]; |
| 1449 | if (iRow >= 0) |
| 1450 | newIndices[newSize++] = iRow; |
| 1451 | } |
| 1452 | } |
| 1453 | startPositive_[numberColumns_] = newSize; |
| 1454 | delete[] which; |
| 1455 | delete[] indices_; |
| 1456 | indices_ = newIndices; |
| 1457 | numberRows_ = newNumber; |
| 1458 | } |
| 1459 | bool ClpPlusMinusOneMatrix::isColOrdered() const |
| 1460 | { |
| 1461 | return columnOrdered_; |
1532 | | if (startNegative_[i] < startPositive_[i+1]) |
1533 | | otherFlags_ &= ~1; // not all +1 |
1534 | | #endif |
1535 | | if(startNegative_[i] < last) |
1536 | | bad++; |
1537 | | else |
1538 | | last = startNegative_[i]; |
1539 | | } |
1540 | | if(startPositive_[number] < last) |
1541 | | bad++; |
1542 | | CoinAssertHint(!bad, "starts are not monotonic"); |
1543 | | for (CoinBigIndex cbi = 0; cbi < numberElements; cbi++) { |
1544 | | maxIndex = CoinMax(indices_[cbi], maxIndex); |
1545 | | minIndex = CoinMin(indices_[cbi], minIndex); |
1546 | | } |
1547 | | CoinAssert(maxIndex < (columnOrdered_ ? numberRows_ : numberColumns_)); |
1548 | | CoinAssert(minIndex >= 0); |
1549 | | if (detail) { |
1550 | | if (minIndex > 0 || maxIndex + 1 < (columnOrdered_ ? numberRows_ : numberColumns_)) |
1551 | | printf("Not full range of indices - %d to %d\n", minIndex, maxIndex); |
1552 | | } |
| 1509 | if (startNegative_[i] < startPositive_[i + 1]) |
| 1510 | otherFlags_ &= ~1; // not all +1 |
| 1511 | #endif |
| 1512 | if (startNegative_[i] < last) |
| 1513 | bad++; |
| 1514 | else |
| 1515 | last = startNegative_[i]; |
| 1516 | } |
| 1517 | if (startPositive_[number] < last) |
| 1518 | bad++; |
| 1519 | CoinAssertHint(!bad, "starts are not monotonic"); |
| 1520 | for (CoinBigIndex cbi = 0; cbi < numberElements; cbi++) { |
| 1521 | maxIndex = CoinMax(indices_[cbi], maxIndex); |
| 1522 | minIndex = CoinMin(indices_[cbi], minIndex); |
| 1523 | } |
| 1524 | CoinAssert(maxIndex < (columnOrdered_ ? numberRows_ : numberColumns_)); |
| 1525 | CoinAssert(minIndex >= 0); |
| 1526 | if (detail) { |
| 1527 | if (minIndex > 0 || maxIndex + 1 < (columnOrdered_ ? numberRows_ : numberColumns_)) |
| 1528 | printf("Not full range of indices - %d to %d\n", minIndex, maxIndex); |
| 1529 | } |
1559 | | ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const |
1560 | | { |
1561 | | int numberRows = model->numberRows(); |
1562 | | int numberColumns = model->numberColumns(); |
1563 | | int number = numberRows + numberColumns; |
1564 | | CoinBigIndex * weights = new CoinBigIndex[number]; |
1565 | | int i; |
1566 | | for (i = 0; i < numberColumns; i++) { |
1567 | | CoinBigIndex j; |
1568 | | CoinBigIndex count = 0; |
1569 | | for (j = startPositive_[i]; j < startPositive_[i+1]; j++) { |
1570 | | int iRow = indices_[j]; |
1571 | | count += inputWeights[iRow]; |
1572 | | } |
1573 | | weights[i] = count; |
1574 | | } |
1575 | | for (i = 0; i < numberRows; i++) { |
1576 | | weights[i+numberColumns] = inputWeights[i]; |
1577 | | } |
1578 | | return weights; |
| 1536 | ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex *model, int *inputWeights) const |
| 1537 | { |
| 1538 | int numberRows = model->numberRows(); |
| 1539 | int numberColumns = model->numberColumns(); |
| 1540 | int number = numberRows + numberColumns; |
| 1541 | CoinBigIndex *weights = new CoinBigIndex[number]; |
| 1542 | int i; |
| 1543 | for (i = 0; i < numberColumns; i++) { |
| 1544 | CoinBigIndex j; |
| 1545 | CoinBigIndex count = 0; |
| 1546 | for (j = startPositive_[i]; j < startPositive_[i + 1]; j++) { |
| 1547 | int iRow = indices_[j]; |
| 1548 | count += inputWeights[iRow]; |
| 1549 | } |
| 1550 | weights[i] = count; |
| 1551 | } |
| 1552 | for (i = 0; i < numberRows; i++) { |
| 1553 | weights[i + numberColumns] = inputWeights[i]; |
| 1554 | } |
| 1555 | return weights; |
1589 | | otherFlags_ = 0; |
1590 | | #endif |
1591 | | for (iColumn = 0; iColumn < number; iColumn++) { |
1592 | | int n = columns[iColumn]->getNumElements(); |
1593 | | const double * element = columns[iColumn]->getElements(); |
1594 | | size += n; |
1595 | | int i; |
1596 | | for (i = 0; i < n; i++) { |
1597 | | if (fabs(element[i]) != 1.0) |
1598 | | numberBad++; |
1599 | | } |
1600 | | } |
1601 | | if (numberBad) |
1602 | | throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix"); |
1603 | | // Get rid of temporary arrays |
1604 | | delete [] lengths_; |
1605 | | lengths_ = NULL; |
1606 | | delete matrix_; |
1607 | | matrix_ = NULL; |
1608 | | CoinBigIndex numberNow = startPositive_[numberColumns_]; |
1609 | | CoinBigIndex * temp; |
1610 | | temp = new CoinBigIndex [numberColumns_+1+number]; |
1611 | | CoinMemcpyN(startPositive_, (numberColumns_ + 1), temp); |
1612 | | delete [] startPositive_; |
1613 | | startPositive_ = temp; |
1614 | | temp = new CoinBigIndex [numberColumns_+number]; |
1615 | | CoinMemcpyN(startNegative_, numberColumns_, temp); |
1616 | | delete [] startNegative_; |
1617 | | startNegative_ = temp; |
1618 | | int * temp2 = new int [numberNow+size]; |
1619 | | CoinMemcpyN(indices_, numberNow, temp2); |
1620 | | delete [] indices_; |
1621 | | indices_ = temp2; |
1622 | | // now add |
1623 | | size = numberNow; |
1624 | | for (iColumn = 0; iColumn < number; iColumn++) { |
1625 | | int n = columns[iColumn]->getNumElements(); |
1626 | | const int * row = columns[iColumn]->getIndices(); |
1627 | | const double * element = columns[iColumn]->getElements(); |
1628 | | int i; |
1629 | | for (i = 0; i < n; i++) { |
1630 | | if (element[i] == 1.0) |
1631 | | indices_[size++] = row[i]; |
1632 | | } |
1633 | | startNegative_[iColumn+numberColumns_] = size; |
1634 | | for (i = 0; i < n; i++) { |
1635 | | if (element[i] == -1.0) |
1636 | | indices_[size++] = row[i]; |
1637 | | } |
1638 | | startPositive_[iColumn+numberColumns_+1] = size; |
1639 | | } |
| 1565 | otherFlags_ = 0; |
| 1566 | #endif |
| 1567 | for (iColumn = 0; iColumn < number; iColumn++) { |
| 1568 | int n = columns[iColumn]->getNumElements(); |
| 1569 | const double *element = columns[iColumn]->getElements(); |
| 1570 | size += n; |
| 1571 | int i; |
| 1572 | for (i = 0; i < n; i++) { |
| 1573 | if (fabs(element[i]) != 1.0) |
| 1574 | numberBad++; |
| 1575 | } |
| 1576 | } |
| 1577 | if (numberBad) |
| 1578 | throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix"); |
| 1579 | // Get rid of temporary arrays |
| 1580 | delete[] lengths_; |
| 1581 | lengths_ = NULL; |
| 1582 | delete matrix_; |
| 1583 | matrix_ = NULL; |
| 1584 | CoinBigIndex numberNow = startPositive_[numberColumns_]; |
| 1585 | CoinBigIndex *temp; |
| 1586 | temp = new CoinBigIndex[numberColumns_ + 1 + number]; |
| 1587 | CoinMemcpyN(startPositive_, (numberColumns_ + 1), temp); |
| 1588 | delete[] startPositive_; |
| 1589 | startPositive_ = temp; |
| 1590 | temp = new CoinBigIndex[numberColumns_ + number]; |
| 1591 | CoinMemcpyN(startNegative_, numberColumns_, temp); |
| 1592 | delete[] startNegative_; |
| 1593 | startNegative_ = temp; |
| 1594 | int *temp2 = new int[numberNow + size]; |
| 1595 | CoinMemcpyN(indices_, numberNow, temp2); |
| 1596 | delete[] indices_; |
| 1597 | indices_ = temp2; |
| 1598 | // now add |
| 1599 | size = numberNow; |
| 1600 | for (iColumn = 0; iColumn < number; iColumn++) { |
| 1601 | int n = columns[iColumn]->getNumElements(); |
| 1602 | const int *row = columns[iColumn]->getIndices(); |
| 1603 | const double *element = columns[iColumn]->getElements(); |
| 1604 | int i; |
| 1605 | for (i = 0; i < n; i++) { |
| 1606 | if (element[i] == 1.0) |
| 1607 | indices_[size++] = row[i]; |
| 1608 | } |
| 1609 | startNegative_[iColumn + numberColumns_] = size; |
| 1610 | for (i = 0; i < n; i++) { |
| 1611 | if (element[i] == -1.0) |
| 1612 | indices_[size++] = row[i]; |
| 1613 | } |
| 1614 | startPositive_[iColumn + numberColumns_ + 1] = size; |
| 1615 | } |
1644 | | void |
1645 | | ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows) |
1646 | | { |
1647 | | // Allocate arrays to use for counting |
1648 | | int * countPositive = new int [numberColumns_+1]; |
1649 | | memset(countPositive, 0, numberColumns_ * sizeof(int)); |
1650 | | int * countNegative = new int [numberColumns_]; |
1651 | | memset(countNegative, 0, numberColumns_ * sizeof(int)); |
1652 | | int iRow; |
1653 | | CoinBigIndex size = 0; |
1654 | | int numberBad = 0; |
1655 | | // say if all +1 |
| 1620 | void ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase *const *rows) |
| 1621 | { |
| 1622 | // Allocate arrays to use for counting |
| 1623 | int *countPositive = new int[numberColumns_ + 1]; |
| 1624 | memset(countPositive, 0, numberColumns_ * sizeof(int)); |
| 1625 | int *countNegative = new int[numberColumns_]; |
| 1626 | memset(countNegative, 0, numberColumns_ * sizeof(int)); |
| 1627 | int iRow; |
| 1628 | CoinBigIndex size = 0; |
| 1629 | int numberBad = 0; |
| 1630 | // say if all +1 |
1657 | | otherFlags_ = 0; |
1658 | | #endif |
1659 | | for (iRow = 0; iRow < number; iRow++) { |
1660 | | int n = rows[iRow]->getNumElements(); |
1661 | | const int * column = rows[iRow]->getIndices(); |
1662 | | const double * element = rows[iRow]->getElements(); |
1663 | | size += n; |
1664 | | int i; |
1665 | | for (i = 0; i < n; i++) { |
1666 | | int iColumn = column[i]; |
1667 | | if (element[i] == 1.0) |
1668 | | countPositive[iColumn]++; |
1669 | | else if (element[i] == -1.0) |
1670 | | countNegative[iColumn]++; |
1671 | | else |
1672 | | numberBad++; |
1673 | | } |
1674 | | } |
1675 | | if (numberBad) |
1676 | | throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix"); |
1677 | | // Get rid of temporary arrays |
1678 | | delete [] lengths_; |
1679 | | lengths_ = NULL; |
1680 | | delete matrix_; |
1681 | | matrix_ = NULL; |
1682 | | CoinBigIndex numberNow = startPositive_[numberColumns_]; |
1683 | | int * newIndices = new int [numberNow+size]; |
1684 | | // Update starts and turn counts into positions |
1685 | | // also move current indices |
1686 | | int iColumn; |
1687 | | CoinBigIndex numberAdded = 0; |
1688 | | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
1689 | | int n, move; |
1690 | | CoinBigIndex now; |
1691 | | now = startPositive_[iColumn]; |
1692 | | move = static_cast<int>(startNegative_[iColumn] - now); |
1693 | | n = countPositive[iColumn]; |
1694 | | startPositive_[iColumn] += numberAdded; |
1695 | | CoinMemcpyN(indices_+now, move, newIndices + startPositive_[iColumn]); |
1696 | | countPositive[iColumn] = static_cast<int>(startNegative_[iColumn] + numberAdded); |
1697 | | numberAdded += n; |
1698 | | now = startNegative_[iColumn]; |
1699 | | move = static_cast<int>(startPositive_[iColumn+1] - now); |
1700 | | n = countNegative[iColumn]; |
1701 | | startNegative_[iColumn] += numberAdded; |
1702 | | CoinMemcpyN(indices_+now, move, newIndices + startNegative_[iColumn]); |
1703 | | countNegative[iColumn] = static_cast<int>(startPositive_[iColumn+1] + numberAdded); |
1704 | | numberAdded += n; |
1705 | | } |
1706 | | delete [] indices_; |
1707 | | indices_ = newIndices; |
1708 | | startPositive_[numberColumns_] += numberAdded; |
1709 | | // Now put in |
1710 | | for (iRow = 0; iRow < number; iRow++) { |
1711 | | int newRow = numberRows_ + iRow; |
1712 | | int n = rows[iRow]->getNumElements(); |
1713 | | const int * column = rows[iRow]->getIndices(); |
1714 | | const double * element = rows[iRow]->getElements(); |
1715 | | int i; |
1716 | | for (i = 0; i < n; i++) { |
1717 | | int iColumn = column[i]; |
1718 | | int put; |
1719 | | if (element[i] == 1.0) { |
1720 | | put = countPositive[iColumn]; |
1721 | | countPositive[iColumn] = put + 1; |
1722 | | } else { |
1723 | | put = countNegative[iColumn]; |
1724 | | countNegative[iColumn] = put + 1; |
1725 | | } |
1726 | | indices_[put] = newRow; |
1727 | | } |
1728 | | } |
1729 | | delete [] countPositive; |
1730 | | delete [] countNegative; |
1731 | | numberRows_ += number; |
| 1632 | otherFlags_ = 0; |
| 1633 | #endif |
| 1634 | for (iRow = 0; iRow < number; iRow++) { |
| 1635 | int n = rows[iRow]->getNumElements(); |
| 1636 | const int *column = rows[iRow]->getIndices(); |
| 1637 | const double *element = rows[iRow]->getElements(); |
| 1638 | size += n; |
| 1639 | int i; |
| 1640 | for (i = 0; i < n; i++) { |
| 1641 | int iColumn = column[i]; |
| 1642 | if (element[i] == 1.0) |
| 1643 | countPositive[iColumn]++; |
| 1644 | else if (element[i] == -1.0) |
| 1645 | countNegative[iColumn]++; |
| 1646 | else |
| 1647 | numberBad++; |
| 1648 | } |
| 1649 | } |
| 1650 | if (numberBad) |
| 1651 | throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix"); |
| 1652 | // Get rid of temporary arrays |
| 1653 | delete[] lengths_; |
| 1654 | lengths_ = NULL; |
| 1655 | delete matrix_; |
| 1656 | matrix_ = NULL; |
| 1657 | CoinBigIndex numberNow = startPositive_[numberColumns_]; |
| 1658 | int *newIndices = new int[numberNow + size]; |
| 1659 | // Update starts and turn counts into positions |
| 1660 | // also move current indices |
| 1661 | int iColumn; |
| 1662 | CoinBigIndex numberAdded = 0; |
| 1663 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
| 1664 | int n, move; |
| 1665 | CoinBigIndex now; |
| 1666 | now = startPositive_[iColumn]; |
| 1667 | move = static_cast< int >(startNegative_[iColumn] - now); |
| 1668 | n = countPositive[iColumn]; |
| 1669 | startPositive_[iColumn] += numberAdded; |
| 1670 | CoinMemcpyN(indices_ + now, move, newIndices + startPositive_[iColumn]); |
| 1671 | countPositive[iColumn] = static_cast< int >(startNegative_[iColumn] + numberAdded); |
| 1672 | numberAdded += n; |
| 1673 | now = startNegative_[iColumn]; |
| 1674 | move = static_cast< int >(startPositive_[iColumn + 1] - now); |
| 1675 | n = countNegative[iColumn]; |
| 1676 | startNegative_[iColumn] += numberAdded; |
| 1677 | CoinMemcpyN(indices_ + now, move, newIndices + startNegative_[iColumn]); |
| 1678 | countNegative[iColumn] = static_cast< int >(startPositive_[iColumn + 1] + numberAdded); |
| 1679 | numberAdded += n; |
| 1680 | } |
| 1681 | delete[] indices_; |
| 1682 | indices_ = newIndices; |
| 1683 | startPositive_[numberColumns_] += numberAdded; |
| 1684 | // Now put in |
| 1685 | for (iRow = 0; iRow < number; iRow++) { |
| 1686 | int newRow = numberRows_ + iRow; |
| 1687 | int n = rows[iRow]->getNumElements(); |
| 1688 | const int *column = rows[iRow]->getIndices(); |
| 1689 | const double *element = rows[iRow]->getElements(); |
| 1690 | int i; |
| 1691 | for (i = 0; i < n; i++) { |
| 1692 | int iColumn = column[i]; |
| 1693 | int put; |
| 1694 | if (element[i] == 1.0) { |
| 1695 | put = countPositive[iColumn]; |
| 1696 | countPositive[iColumn] = put + 1; |
| 1697 | } else { |
| 1698 | put = countNegative[iColumn]; |
| 1699 | countNegative[iColumn] = put + 1; |
| 1700 | } |
| 1701 | indices_[put] = newRow; |
| 1702 | } |
| 1703 | } |
| 1704 | delete[] countPositive; |
| 1705 | delete[] countNegative; |
| 1706 | numberRows_ += number; |
1771 | | void |
1772 | | ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, |
1773 | | int & bestSequence, int & numberWanted) |
1774 | | { |
1775 | | numberWanted = currentWanted_; |
1776 | | int start = static_cast<int> (startFraction * numberColumns_); |
1777 | | int end = CoinMin(static_cast<int> (endFraction * numberColumns_ + 1), numberColumns_); |
1778 | | CoinBigIndex j; |
1779 | | double tolerance = model->currentDualTolerance(); |
1780 | | double * COIN_RESTRICT reducedCost = model->djRegion(); |
1781 | | const double * COIN_RESTRICT duals = model->dualRowSolution(); |
1782 | | const double * COIN_RESTRICT cost = model->costRegion(); |
1783 | | double bestDj; |
1784 | | if (bestSequence >= 0) |
1785 | | bestDj = fabs(reducedCost[bestSequence]); |
1786 | | else |
1787 | | bestDj = tolerance; |
1788 | | int sequenceOut = model->sequenceOut(); |
1789 | | int saveSequence = bestSequence; |
1790 | | int iSequence; |
1791 | | for (iSequence = start; iSequence < end; iSequence++) { |
1792 | | if (iSequence != sequenceOut) { |
1793 | | double value; |
1794 | | ClpSimplex::Status status = model->getStatus(iSequence); |
| 1744 | void ClpPlusMinusOneMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction, |
| 1745 | int &bestSequence, int &numberWanted) |
| 1746 | { |
| 1747 | numberWanted = currentWanted_; |
| 1748 | int start = static_cast< int >(startFraction * numberColumns_); |
| 1749 | int end = CoinMin(static_cast< int >(endFraction * numberColumns_ + 1), numberColumns_); |
| 1750 | CoinBigIndex j; |
| 1751 | double tolerance = model->currentDualTolerance(); |
| 1752 | double *COIN_RESTRICT reducedCost = model->djRegion(); |
| 1753 | const double *COIN_RESTRICT duals = model->dualRowSolution(); |
| 1754 | const double *COIN_RESTRICT cost = model->costRegion(); |
| 1755 | double bestDj; |
| 1756 | if (bestSequence >= 0) |
| 1757 | bestDj = fabs(reducedCost[bestSequence]); |
| 1758 | else |
| 1759 | bestDj = tolerance; |
| 1760 | int sequenceOut = model->sequenceOut(); |
| 1761 | int saveSequence = bestSequence; |
| 1762 | int iSequence; |
| 1763 | for (iSequence = start; iSequence < end; iSequence++) { |
| 1764 | if (iSequence != sequenceOut) { |
| 1765 | double value; |
| 1766 | ClpSimplex::Status status = model->getStatus(iSequence); |
1798 | | case ClpSimplex::basic: |
1799 | | case ClpSimplex::isFixed: |
1800 | | break; |
1801 | | case ClpSimplex::isFree: |
1802 | | case ClpSimplex::superBasic: |
1803 | | value = cost[iSequence]; |
1804 | | j = startPositive_[iSequence]; |
1805 | | for (; j < startNegative_[iSequence]; j++) { |
1806 | | int iRow = indices_[j]; |
1807 | | value -= duals[iRow]; |
1808 | | } |
1809 | | for (; j < startPositive_[iSequence+1]; j++) { |
1810 | | int iRow = indices_[j]; |
1811 | | value += duals[iRow]; |
1812 | | } |
1813 | | value = fabs(value); |
1814 | | if (value > FREE_ACCEPT * tolerance) { |
1815 | | numberWanted--; |
1816 | | // we are going to bias towards free (but only if reasonable) |
1817 | | value *= FREE_BIAS; |
1818 | | if (value > bestDj) { |
1819 | | // check flagged variable and correct dj |
1820 | | if (!model->flagged(iSequence)) { |
1821 | | bestDj = value; |
1822 | | bestSequence = iSequence; |
1823 | | } else { |
1824 | | // just to make sure we don't exit before got something |
1825 | | numberWanted++; |
1826 | | } |
1827 | | } |
1828 | | } |
1829 | | break; |
1830 | | case ClpSimplex::atUpperBound: |
1831 | | value = cost[iSequence]; |
1832 | | j = startPositive_[iSequence]; |
1833 | | for (; j < startNegative_[iSequence]; j++) { |
1834 | | int iRow = indices_[j]; |
1835 | | value -= duals[iRow]; |
1836 | | } |
1837 | | for (; j < startPositive_[iSequence+1]; j++) { |
1838 | | int iRow = indices_[j]; |
1839 | | value += duals[iRow]; |
1840 | | } |
1841 | | if (value > tolerance) { |
1842 | | numberWanted--; |
1843 | | if (value > bestDj) { |
1844 | | // check flagged variable and correct dj |
1845 | | if (!model->flagged(iSequence)) { |
1846 | | bestDj = value; |
1847 | | bestSequence = iSequence; |
1848 | | } else { |
1849 | | // just to make sure we don't exit before got something |
1850 | | numberWanted++; |
1851 | | } |
1852 | | } |
1853 | | } |
1854 | | break; |
1855 | | case ClpSimplex::atLowerBound: |
1856 | | value = cost[iSequence]; |
1857 | | j = startPositive_[iSequence]; |
1858 | | for (; j < startNegative_[iSequence]; j++) { |
1859 | | int iRow = indices_[j]; |
1860 | | value -= duals[iRow]; |
1861 | | } |
1862 | | for (; j < startPositive_[iSequence+1]; j++) { |
1863 | | int iRow = indices_[j]; |
1864 | | value += duals[iRow]; |
1865 | | } |
1866 | | value = -value; |
1867 | | if (value > tolerance) { |
1868 | | numberWanted--; |
1869 | | if (value > bestDj) { |
1870 | | // check flagged variable and correct dj |
1871 | | if (!model->flagged(iSequence)) { |
1872 | | bestDj = value; |
1873 | | bestSequence = iSequence; |
1874 | | } else { |
1875 | | // just to make sure we don't exit before got something |
1876 | | numberWanted++; |
1877 | | } |
1878 | | } |
1879 | | } |
1880 | | break; |
1881 | | } |
| 1770 | case ClpSimplex::basic: |
| 1771 | case ClpSimplex::isFixed: |
| 1772 | break; |
| 1773 | case ClpSimplex::isFree: |
| 1774 | case ClpSimplex::superBasic: |
| 1775 | value = cost[iSequence]; |
| 1776 | j = startPositive_[iSequence]; |
| 1777 | for (; j < startNegative_[iSequence]; j++) { |
| 1778 | int iRow = indices_[j]; |
| 1779 | value -= duals[iRow]; |
| 1780 | } |
| 1781 | for (; j < startPositive_[iSequence + 1]; j++) { |
| 1782 | int iRow = indices_[j]; |
| 1783 | value += duals[iRow]; |
| 1784 | } |
| 1785 | value = fabs(value); |
| 1786 | if (value > FREE_ACCEPT * tolerance) { |
| 1787 | numberWanted--; |
| 1788 | // we are going to bias towards free (but only if reasonable) |
| 1789 | value *= FREE_BIAS; |
| 1790 | if (value > bestDj) { |
| 1791 | // check flagged variable and correct dj |
| 1792 | if (!model->flagged(iSequence)) { |
| 1793 | bestDj = value; |
| 1794 | bestSequence = iSequence; |
| 1795 | } else { |
| 1796 | // just to make sure we don't exit before got something |
| 1797 | numberWanted++; |
| 1798 | } |
1919 | | bool |
1920 | | ClpPlusMinusOneMatrix::canCombine(const ClpSimplex * model, |
1921 | | const CoinIndexedVector * pi) const |
1922 | | { |
1923 | | int numberInRowArray = pi->getNumElements(); |
1924 | | int numberRows = model->numberRows(); |
1925 | | bool packed = pi->packedMode(); |
1926 | | // factor should be smaller if doing both with two pi vectors |
1927 | | double factor = 0.27; |
1928 | | // We may not want to do by row if there may be cache problems |
1929 | | // It would be nice to find L2 cache size - for moment 512K |
1930 | | // Be slightly optimistic |
1931 | | if (numberColumns_ * sizeof(double) > 1000000) { |
1932 | | if (numberRows * 10 < numberColumns_) |
1933 | | factor *= 0.333333333; |
1934 | | else if (numberRows * 4 < numberColumns_) |
1935 | | factor *= 0.5; |
1936 | | else if (numberRows * 2 < numberColumns_) |
1937 | | factor *= 0.66666666667; |
1938 | | //if (model->numberIterations()%50==0) |
1939 | | //printf("%d nonzero\n",numberInRowArray); |
1940 | | } |
1941 | | // if not packed then bias a bit more towards by column |
1942 | | if (!packed) |
1943 | | factor *= 0.9; |
1944 | | return (numberInRowArray > factor * numberRows || !model->rowCopy()); |
| 1890 | bool ClpPlusMinusOneMatrix::canCombine(const ClpSimplex *model, |
| 1891 | const CoinIndexedVector *pi) const |
| 1892 | { |
| 1893 | int numberInRowArray = pi->getNumElements(); |
| 1894 | int numberRows = model->numberRows(); |
| 1895 | bool packed = pi->packedMode(); |
| 1896 | // factor should be smaller if doing both with two pi vectors |
| 1897 | double factor = 0.27; |
| 1898 | // We may not want to do by row if there may be cache problems |
| 1899 | // It would be nice to find L2 cache size - for moment 512K |
| 1900 | // Be slightly optimistic |
| 1901 | if (numberColumns_ * sizeof(double) > 1000000) { |
| 1902 | if (numberRows * 10 < numberColumns_) |
| 1903 | factor *= 0.333333333; |
| 1904 | else if (numberRows * 4 < numberColumns_) |
| 1905 | factor *= 0.5; |
| 1906 | else if (numberRows * 2 < numberColumns_) |
| 1907 | factor *= 0.66666666667; |
| 1908 | //if (model->numberIterations()%50==0) |
| 1909 | //printf("%d nonzero\n",numberInRowArray); |
| 1910 | } |
| 1911 | // if not packed then bias a bit more towards by column |
| 1912 | if (!packed) |
| 1913 | factor *= 0.9; |
| 1914 | return (numberInRowArray > factor * numberRows || !model->rowCopy()); |
1951 | | int |
1952 | | ClpPlusMinusOneMatrix::transposeTimes2(const ClpSimplex * model, |
1953 | | const CoinIndexedVector * pi1, CoinIndexedVector * dj1, |
1954 | | const CoinIndexedVector * pi2, |
1955 | | CoinIndexedVector * spare, |
1956 | | double * , double * , |
1957 | | double referenceIn, double devex, |
1958 | | // Array for exact devex to say what is in reference framework |
1959 | | unsigned int * COIN_RESTRICT reference, |
1960 | | double * COIN_RESTRICT weights, double scaleFactor) |
1961 | | { |
1962 | | // put row of tableau in dj1 |
1963 | | double * COIN_RESTRICT pi = pi1->denseVector(); |
1964 | | int numberNonZero = 0; |
1965 | | int * COIN_RESTRICT index = dj1->getIndices(); |
1966 | | double * COIN_RESTRICT array = dj1->denseVector(); |
1967 | | int numberInRowArray = pi1->getNumElements(); |
1968 | | double zeroTolerance = model->zeroTolerance(); |
1969 | | bool packed = pi1->packedMode(); |
1970 | | // do by column |
1971 | | int iColumn; |
1972 | | assert (!spare->getNumElements()); |
1973 | | double * COIN_RESTRICT piWeight = pi2->denseVector(); |
1974 | | assert (!pi2->packedMode()); |
1975 | | bool killDjs = (scaleFactor == 0.0); |
1976 | | if (!scaleFactor) |
1977 | | scaleFactor = 1.0; |
1978 | | // Note scale factor was -1.0 |
1979 | | if (packed) { |
1980 | | // need to expand pi into y |
1981 | | assert(spare->capacity() >= model->numberRows()); |
1982 | | double * COIN_RESTRICT piOld = pi; |
1983 | | pi = spare->denseVector(); |
1984 | | const int * COIN_RESTRICT whichRow = pi1->getIndices(); |
1985 | | int i; |
1986 | | // modify pi so can collapse to one loop |
1987 | | for (i = 0; i < numberInRowArray; i++) { |
1988 | | int iRow = whichRow[i]; |
1989 | | pi[iRow] = piOld[i]; |
| 1921 | int ClpPlusMinusOneMatrix::transposeTimes2(const ClpSimplex *model, |
| 1922 | const CoinIndexedVector *pi1, CoinIndexedVector *dj1, |
| 1923 | const CoinIndexedVector *pi2, |
| 1924 | CoinIndexedVector *spare, |
| 1925 | double *, double *, |
| 1926 | double referenceIn, double devex, |
| 1927 | // Array for exact devex to say what is in reference framework |
| 1928 | unsigned int *COIN_RESTRICT reference, |
| 1929 | double *COIN_RESTRICT weights, double scaleFactor) |
| 1930 | { |
| 1931 | // put row of tableau in dj1 |
| 1932 | double *COIN_RESTRICT pi = pi1->denseVector(); |
| 1933 | int numberNonZero = 0; |
| 1934 | int *COIN_RESTRICT index = dj1->getIndices(); |
| 1935 | double *COIN_RESTRICT array = dj1->denseVector(); |
| 1936 | int numberInRowArray = pi1->getNumElements(); |
| 1937 | double zeroTolerance = model->zeroTolerance(); |
| 1938 | bool packed = pi1->packedMode(); |
| 1939 | // do by column |
| 1940 | int iColumn; |
| 1941 | assert(!spare->getNumElements()); |
| 1942 | double *COIN_RESTRICT piWeight = pi2->denseVector(); |
| 1943 | assert(!pi2->packedMode()); |
| 1944 | bool killDjs = (scaleFactor == 0.0); |
| 1945 | if (!scaleFactor) |
| 1946 | scaleFactor = 1.0; |
| 1947 | // Note scale factor was -1.0 |
| 1948 | if (packed) { |
| 1949 | // need to expand pi into y |
| 1950 | assert(spare->capacity() >= model->numberRows()); |
| 1951 | double *COIN_RESTRICT piOld = pi; |
| 1952 | pi = spare->denseVector(); |
| 1953 | const int *COIN_RESTRICT whichRow = pi1->getIndices(); |
| 1954 | int i; |
| 1955 | // modify pi so can collapse to one loop |
| 1956 | for (i = 0; i < numberInRowArray; i++) { |
| 1957 | int iRow = whichRow[i]; |
| 1958 | pi[iRow] = piOld[i]; |
| 1959 | } |
| 1960 | CoinBigIndex j; |
| 1961 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
| 1962 | ClpSimplex::Status status = model->getStatus(iColumn); |
| 1963 | if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) |
| 1964 | continue; |
| 1965 | double value = 0.0; |
| 1966 | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
| 1967 | int iRow = indices_[j]; |
| 1968 | value -= pi[iRow]; |
| 1969 | } |
| 1970 | for (; j < startPositive_[iColumn + 1]; j++) { |
| 1971 | int iRow = indices_[j]; |
| 1972 | value += pi[iRow]; |
| 1973 | } |
| 1974 | if (fabs(value) > zeroTolerance) { |
| 1975 | // and do other array |
| 1976 | double modification = 0.0; |
| 1977 | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
| 1978 | int iRow = indices_[j]; |
| 1979 | modification += piWeight[iRow]; |
| 1980 | } |
| 1981 | for (; j < startPositive_[iColumn + 1]; j++) { |
| 1982 | int iRow = indices_[j]; |
| 1983 | modification -= piWeight[iRow]; |
| 1984 | } |
| 1985 | double thisWeight = weights[iColumn]; |
| 1986 | double pivot = value * scaleFactor; |
| 1987 | double pivotSquared = pivot * pivot; |
| 1988 | thisWeight += pivotSquared * devex + pivot * modification; |
| 1989 | if (thisWeight < DEVEX_TRY_NORM) { |
| 1990 | if (referenceIn < 0.0) { |
| 1991 | // steepest |
| 1992 | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
| 1993 | } else { |
| 1994 | // exact |
| 1995 | thisWeight = referenceIn * pivotSquared; |
| 1996 | if (reference(iColumn)) |
| 1997 | thisWeight += 1.0; |
| 1998 | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
1991 | | CoinBigIndex j; |
1992 | | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
1993 | | ClpSimplex::Status status = model->getStatus(iColumn); |
1994 | | if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue; |
1995 | | double value = 0.0; |
1996 | | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
1997 | | int iRow = indices_[j]; |
1998 | | value -= pi[iRow]; |
1999 | | } |
2000 | | for (; j < startPositive_[iColumn+1]; j++) { |
2001 | | int iRow = indices_[j]; |
2002 | | value += pi[iRow]; |
2003 | | } |
2004 | | if (fabs(value) > zeroTolerance) { |
2005 | | // and do other array |
2006 | | double modification = 0.0; |
2007 | | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
2008 | | int iRow = indices_[j]; |
2009 | | modification += piWeight[iRow]; |
2010 | | } |
2011 | | for (; j < startPositive_[iColumn+1]; j++) { |
2012 | | int iRow = indices_[j]; |
2013 | | modification -= piWeight[iRow]; |
2014 | | } |
2015 | | double thisWeight = weights[iColumn]; |
2016 | | double pivot = value * scaleFactor; |
2017 | | double pivotSquared = pivot * pivot; |
2018 | | thisWeight += pivotSquared * devex + pivot * modification; |
2019 | | if (thisWeight < DEVEX_TRY_NORM) { |
2020 | | if (referenceIn < 0.0) { |
2021 | | // steepest |
2022 | | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
2023 | | } else { |
2024 | | // exact |
2025 | | thisWeight = referenceIn * pivotSquared; |
2026 | | if (reference(iColumn)) |
2027 | | thisWeight += 1.0; |
2028 | | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
2029 | | } |
2030 | | } |
2031 | | weights[iColumn] = thisWeight; |
2032 | | if (!killDjs) { |
2033 | | array[numberNonZero] = value; |
2034 | | index[numberNonZero++] = iColumn; |
2035 | | } |
2036 | | } |
| 2000 | } |
| 2001 | weights[iColumn] = thisWeight; |
| 2002 | if (!killDjs) { |
| 2003 | array[numberNonZero] = value; |
| 2004 | index[numberNonZero++] = iColumn; |
| 2005 | } |
| 2006 | } |
| 2007 | } |
| 2008 | // zero out |
| 2009 | for (i = 0; i < numberInRowArray; i++) { |
| 2010 | int iRow = whichRow[i]; |
| 2011 | pi[iRow] = 0.0; |
| 2012 | } |
| 2013 | } else { |
| 2014 | CoinBigIndex j; |
| 2015 | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
| 2016 | ClpSimplex::Status status = model->getStatus(iColumn); |
| 2017 | if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) |
| 2018 | continue; |
| 2019 | double value = 0.0; |
| 2020 | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
| 2021 | int iRow = indices_[j]; |
| 2022 | value -= pi[iRow]; |
| 2023 | } |
| 2024 | for (; j < startPositive_[iColumn + 1]; j++) { |
| 2025 | int iRow = indices_[j]; |
| 2026 | value += pi[iRow]; |
| 2027 | } |
| 2028 | if (fabs(value) > zeroTolerance) { |
| 2029 | // and do other array |
| 2030 | double modification = 0.0; |
| 2031 | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
| 2032 | int iRow = indices_[j]; |
| 2033 | modification += piWeight[iRow]; |
| 2034 | } |
| 2035 | for (; j < startPositive_[iColumn + 1]; j++) { |
| 2036 | int iRow = indices_[j]; |
| 2037 | modification -= piWeight[iRow]; |
| 2038 | } |
| 2039 | double thisWeight = weights[iColumn]; |
| 2040 | double pivot = value * scaleFactor; |
| 2041 | double pivotSquared = pivot * pivot; |
| 2042 | thisWeight += pivotSquared * devex + pivot * modification; |
| 2043 | if (thisWeight < DEVEX_TRY_NORM) { |
| 2044 | if (referenceIn < 0.0) { |
| 2045 | // steepest |
| 2046 | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
| 2047 | } else { |
| 2048 | // exact |
| 2049 | thisWeight = referenceIn * pivotSquared; |
| 2050 | if (reference(iColumn)) |
| 2051 | thisWeight += 1.0; |
| 2052 | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
2038 | | // zero out |
2039 | | for (i = 0; i < numberInRowArray; i++) { |
2040 | | int iRow = whichRow[i]; |
2041 | | pi[iRow] = 0.0; |
2042 | | } |
2043 | | } else { |
2044 | | CoinBigIndex j; |
2045 | | for (iColumn = 0; iColumn < numberColumns_; iColumn++) { |
2046 | | ClpSimplex::Status status = model->getStatus(iColumn); |
2047 | | if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue; |
2048 | | double value = 0.0; |
2049 | | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
2050 | | int iRow = indices_[j]; |
2051 | | value -= pi[iRow]; |
2052 | | } |
2053 | | for (; j < startPositive_[iColumn+1]; j++) { |
2054 | | int iRow = indices_[j]; |
2055 | | value += pi[iRow]; |
2056 | | } |
2057 | | if (fabs(value) > zeroTolerance) { |
2058 | | // and do other array |
2059 | | double modification = 0.0; |
2060 | | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
2061 | | int iRow = indices_[j]; |
2062 | | modification += piWeight[iRow]; |
2063 | | } |
2064 | | for (; j < startPositive_[iColumn+1]; j++) { |
2065 | | int iRow = indices_[j]; |
2066 | | modification -= piWeight[iRow]; |
2067 | | } |
2068 | | double thisWeight = weights[iColumn]; |
2069 | | double pivot = value * scaleFactor; |
2070 | | double pivotSquared = pivot * pivot; |
2071 | | thisWeight += pivotSquared * devex + pivot * modification; |
2072 | | if (thisWeight < DEVEX_TRY_NORM) { |
2073 | | if (referenceIn < 0.0) { |
2074 | | // steepest |
2075 | | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
2076 | | } else { |
2077 | | // exact |
2078 | | thisWeight = referenceIn * pivotSquared; |
2079 | | if (reference(iColumn)) |
2080 | | thisWeight += 1.0; |
2081 | | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
2082 | | } |
2083 | | } |
2084 | | weights[iColumn] = thisWeight; |
2085 | | if (!killDjs) { |
2086 | | array[iColumn] = value; |
2087 | | index[numberNonZero++] = iColumn; |
2088 | | } |
2089 | | } |
2090 | | } |
2091 | | } |
2092 | | dj1->setNumElements(numberNonZero); |
2093 | | spare->setNumElements(0); |
2094 | | if (packed) |
2095 | | dj1->setPackedMode(true); |
2096 | | return 0; |
| 2054 | } |
| 2055 | weights[iColumn] = thisWeight; |
| 2056 | if (!killDjs) { |
| 2057 | array[iColumn] = value; |
| 2058 | index[numberNonZero++] = iColumn; |
| 2059 | } |
| 2060 | } |
| 2061 | } |
| 2062 | } |
| 2063 | dj1->setNumElements(numberNonZero); |
| 2064 | spare->setNumElements(0); |
| 2065 | if (packed) |
| 2066 | dj1->setPackedMode(true); |
| 2067 | return 0; |
2099 | | void |
2100 | | ClpPlusMinusOneMatrix::subsetTimes2(const ClpSimplex * , |
2101 | | CoinIndexedVector * dj1, |
2102 | | const CoinIndexedVector * pi2, CoinIndexedVector *, |
2103 | | double referenceIn, double devex, |
2104 | | // Array for exact devex to say what is in reference framework |
2105 | | unsigned int * COIN_RESTRICT reference, |
2106 | | double * COIN_RESTRICT weights, double scaleFactor) |
2107 | | { |
2108 | | int number = dj1->getNumElements(); |
2109 | | const int * COIN_RESTRICT index = dj1->getIndices(); |
2110 | | double * COIN_RESTRICT array = dj1->denseVector(); |
2111 | | assert( dj1->packedMode()); |
| 2070 | void ClpPlusMinusOneMatrix::subsetTimes2(const ClpSimplex *, |
| 2071 | CoinIndexedVector *dj1, |
| 2072 | const CoinIndexedVector *pi2, CoinIndexedVector *, |
| 2073 | double referenceIn, double devex, |
| 2074 | // Array for exact devex to say what is in reference framework |
| 2075 | unsigned int *COIN_RESTRICT reference, |
| 2076 | double *COIN_RESTRICT weights, double scaleFactor) |
| 2077 | { |
| 2078 | int number = dj1->getNumElements(); |
| 2079 | const int *COIN_RESTRICT index = dj1->getIndices(); |
| 2080 | double *COIN_RESTRICT array = dj1->denseVector(); |
| 2081 | assert(dj1->packedMode()); |
2113 | | double * COIN_RESTRICT piWeight = pi2->denseVector(); |
2114 | | bool killDjs = (scaleFactor == 0.0); |
2115 | | if (!scaleFactor) |
2116 | | scaleFactor = 1.0; |
2117 | | for (int k = 0; k < number; k++) { |
2118 | | int iColumn = index[k]; |
2119 | | double pivot = array[k] * scaleFactor; |
2120 | | if (killDjs) |
2121 | | array[k] = 0.0; |
2122 | | // and do other array |
2123 | | double modification = 0.0; |
2124 | | CoinBigIndex j; |
2125 | | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
2126 | | int iRow = indices_[j]; |
2127 | | modification += piWeight[iRow]; |
2128 | | } |
2129 | | for (; j < startPositive_[iColumn+1]; j++) { |
2130 | | int iRow = indices_[j]; |
2131 | | modification -= piWeight[iRow]; |
2132 | | } |
2133 | | double thisWeight = weights[iColumn]; |
2134 | | double pivotSquared = pivot * pivot; |
2135 | | thisWeight += pivotSquared * devex + pivot * modification; |
2136 | | if (thisWeight < DEVEX_TRY_NORM) { |
2137 | | if (referenceIn < 0.0) { |
2138 | | // steepest |
2139 | | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
2140 | | } else { |
2141 | | // exact |
2142 | | thisWeight = referenceIn * pivotSquared; |
2143 | | if (reference(iColumn)) |
2144 | | thisWeight += 1.0; |
2145 | | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
2146 | | } |
2147 | | } |
2148 | | weights[iColumn] = thisWeight; |
2149 | | } |
| 2083 | double *COIN_RESTRICT piWeight = pi2->denseVector(); |
| 2084 | bool killDjs = (scaleFactor == 0.0); |
| 2085 | if (!scaleFactor) |
| 2086 | scaleFactor = 1.0; |
| 2087 | for (int k = 0; k < number; k++) { |
| 2088 | int iColumn = index[k]; |
| 2089 | double pivot = array[k] * scaleFactor; |
| 2090 | if (killDjs) |
| 2091 | array[k] = 0.0; |
| 2092 | // and do other array |
| 2093 | double modification = 0.0; |
| 2094 | CoinBigIndex j; |
| 2095 | for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) { |
| 2096 | int iRow = indices_[j]; |
| 2097 | modification += piWeight[iRow]; |
| 2098 | } |
| 2099 | for (; j < startPositive_[iColumn + 1]; j++) { |
| 2100 | int iRow = indices_[j]; |
| 2101 | modification -= piWeight[iRow]; |
| 2102 | } |
| 2103 | double thisWeight = weights[iColumn]; |
| 2104 | double pivotSquared = pivot * pivot; |
| 2105 | thisWeight += pivotSquared * devex + pivot * modification; |
| 2106 | if (thisWeight < DEVEX_TRY_NORM) { |
| 2107 | if (referenceIn < 0.0) { |
| 2108 | // steepest |
| 2109 | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
| 2110 | } else { |
| 2111 | // exact |
| 2112 | thisWeight = referenceIn * pivotSquared; |
| 2113 | if (reference(iColumn)) |
| 2114 | thisWeight += 1.0; |
| 2115 | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
| 2116 | } |
| 2117 | } |
| 2118 | weights[iColumn] = thisWeight; |
| 2119 | } |
2178 | | } else { |
2179 | | length = numberRows_; |
2180 | | numberRows_ = newnumrows; |
2181 | | number = numberRows_; |
2182 | | } |
2183 | | if (number > length) { |
2184 | | CoinBigIndex * temp; |
2185 | | int i; |
2186 | | CoinBigIndex end = startPositive_[length]; |
2187 | | temp = new CoinBigIndex [number+1]; |
2188 | | CoinMemcpyN(startPositive_, (length + 1), temp); |
2189 | | delete [] startPositive_; |
2190 | | for (i = length + 1; i < number + 1; i++) |
2191 | | temp[i] = end; |
2192 | | startPositive_ = temp; |
2193 | | temp = new CoinBigIndex [number]; |
2194 | | CoinMemcpyN(startNegative_, length, temp); |
2195 | | delete [] startNegative_; |
2196 | | for (i = length; i < number; i++) |
2197 | | temp[i] = end; |
2198 | | startNegative_ = temp; |
2199 | | } |
| 2147 | } else { |
| 2148 | length = numberRows_; |
| 2149 | numberRows_ = newnumrows; |
| 2150 | number = numberRows_; |
| 2151 | } |
| 2152 | if (number > length) { |
| 2153 | CoinBigIndex *temp; |
| 2154 | int i; |
| 2155 | CoinBigIndex end = startPositive_[length]; |
| 2156 | temp = new CoinBigIndex[number + 1]; |
| 2157 | CoinMemcpyN(startPositive_, (length + 1), temp); |
| 2158 | delete[] startPositive_; |
| 2159 | for (i = length + 1; i < number + 1; i++) |
| 2160 | temp[i] = end; |
| 2161 | startPositive_ = temp; |
| 2162 | temp = new CoinBigIndex[number]; |
| 2163 | CoinMemcpyN(startNegative_, length, temp); |
| 2164 | delete[] startNegative_; |
| 2165 | for (i = length; i < number; i++) |
| 2166 | temp[i] = end; |
| 2167 | startNegative_ = temp; |
| 2168 | } |
2214 | | otherFlags_ = 0; |
2215 | | #endif |
2216 | | // make into CoinPackedVector |
2217 | | CoinPackedVectorBase ** vectors = |
2218 | | new CoinPackedVectorBase * [number]; |
2219 | | int iVector; |
2220 | | for (iVector = 0; iVector < number; iVector++) { |
2221 | | CoinBigIndex iStart = starts[iVector]; |
2222 | | vectors[iVector] = |
2223 | | new CoinPackedVector(static_cast<int>(starts[iVector+1] - iStart), |
2224 | | index + iStart, element + iStart); |
2225 | | } |
2226 | | if (type == 0) { |
2227 | | // rows |
2228 | | appendRows(number, vectors); |
2229 | | } else { |
2230 | | // columns |
2231 | | appendCols(number, vectors); |
2232 | | } |
2233 | | for (iVector = 0; iVector < number; iVector++) |
2234 | | delete vectors[iVector]; |
2235 | | delete [] vectors; |
2236 | | return numberErrors; |
| 2182 | otherFlags_ = 0; |
| 2183 | #endif |
| 2184 | // make into CoinPackedVector |
| 2185 | CoinPackedVectorBase **vectors = new CoinPackedVectorBase *[number]; |
| 2186 | int iVector; |
| 2187 | for (iVector = 0; iVector < number; iVector++) { |
| 2188 | CoinBigIndex iStart = starts[iVector]; |
| 2189 | vectors[iVector] = new CoinPackedVector(static_cast< int >(starts[iVector + 1] - iStart), |
| 2190 | index + iStart, element + iStart); |
| 2191 | } |
| 2192 | if (type == 0) { |
| 2193 | // rows |
| 2194 | appendRows(number, vectors); |
| 2195 | } else { |
| 2196 | // columns |
| 2197 | appendCols(number, vectors); |
| 2198 | } |
| 2199 | for (iVector = 0; iVector < number; iVector++) |
| 2200 | delete vectors[iVector]; |
| 2201 | delete[] vectors; |
| 2202 | return numberErrors; |
2348 | | while ( true ) { |
2349 | | j1 = hashThis[ipos].index; |
2350 | | if (j1 == -1) { |
2351 | | hashThis[ipos].index=numberDifferent_; |
2352 | | j1=numberDifferent_; |
2353 | | break; |
2354 | | } else if (value==tempDifferent[j1]) { |
2355 | | break; |
2356 | | } else { |
2357 | | int k = hashThis[ipos].next; |
2358 | | if ( k == -1 ) { |
2359 | | j1 = numberDifferent_; |
2360 | | while ( true ) { |
2361 | | ++hashDifferent; |
2362 | | if ( hashThis[hashDifferent].index == -1 ) { |
2363 | | break; |
2364 | | } |
2365 | | } |
2366 | | hashThis[ipos].next = hashDifferent; |
2367 | | hashThis[hashDifferent].index = numberDifferent_; |
2368 | | break; |
2369 | | } else { |
2370 | | ipos = k; |
2371 | | /* nothing worked - try it again */ |
2372 | | } |
2373 | | } |
2374 | | } |
2375 | | #if HASH ==2 |
2376 | | int j=j1; |
2377 | | #endif |
2378 | | #endif |
2379 | | #if HASH ==3 |
2380 | | assert (j==j1); |
2381 | | #endif |
2382 | | if (j==numberDifferent_) { |
2383 | | if (j==maxPool) |
2384 | | break; // too many |
2385 | | tempDifferent[j]=value; |
2386 | | numberDifferent_++; |
2387 | | } |
2388 | | stuff.pool_=j; |
2389 | | stuff_[numberElements++]=stuff; |
2390 | | } |
2391 | | } |
2392 | | columnStart_[numberColumns_]=numberElements; |
2393 | | #if HASH>1 |
2394 | | delete [] hashThis; |
2395 | | #endif |
2396 | | elements_ = new double [numberDifferent_]; |
2397 | | memcpy(elements_,tempDifferent,numberDifferent_*sizeof(double)); |
2398 | | delete [] tempDifferent; |
2399 | | if (numberDifferent_==maxPool) { |
2400 | | delete [] stuff_; |
2401 | | delete [] elements_; |
2402 | | delete [] columnStart_; |
| 2315 | while (true) { |
| 2316 | j1 = hashThis[ipos].index; |
| 2317 | if (j1 == -1) { |
| 2318 | hashThis[ipos].index = numberDifferent_; |
| 2319 | j1 = numberDifferent_; |
| 2320 | break; |
| 2321 | } else if (value == tempDifferent[j1]) { |
| 2322 | break; |
| 2323 | } else { |
| 2324 | int k = hashThis[ipos].next; |
| 2325 | if (k == -1) { |
| 2326 | j1 = numberDifferent_; |
| 2327 | while (true) { |
| 2328 | ++hashDifferent; |
| 2329 | if (hashThis[hashDifferent].index == -1) { |
| 2330 | break; |
| 2331 | } |
| 2332 | } |
| 2333 | hashThis[ipos].next = hashDifferent; |
| 2334 | hashThis[hashDifferent].index = numberDifferent_; |
| 2335 | break; |
| 2336 | } else { |
| 2337 | ipos = k; |
| 2338 | /* nothing worked - try it again */ |
| 2339 | } |
| 2340 | } |
| 2341 | } |
| 2342 | #if HASH == 2 |
| 2343 | int j = j1; |
| 2344 | #endif |
| 2345 | #endif |
| 2346 | #if HASH == 3 |
| 2347 | assert(j == j1); |
| 2348 | #endif |
| 2349 | if (j == numberDifferent_) { |
| 2350 | if (j == maxPool) |
| 2351 | break; // too many |
| 2352 | tempDifferent[j] = value; |
| 2353 | numberDifferent_++; |
| 2354 | } |
| 2355 | stuff.pool_ = j; |
| 2356 | stuff_[numberElements++] = stuff; |
| 2357 | } |
| 2358 | } |
| 2359 | columnStart_[numberColumns_] = numberElements; |
| 2360 | #if HASH > 1 |
| 2361 | delete[] hashThis; |
| 2362 | #endif |
| 2363 | elements_ = new double[numberDifferent_]; |
| 2364 | memcpy(elements_, tempDifferent, numberDifferent_ * sizeof(double)); |
| 2365 | delete[] tempDifferent; |
| 2366 | if (numberDifferent_ == maxPool) { |
| 2367 | delete[] stuff_; |
| 2368 | delete[] elements_; |
| 2369 | delete[] columnStart_; |
2702 | | void |
2703 | | ClpPoolMatrix::times(double scalar, |
2704 | | const double * COIN_RESTRICT x, double * COIN_RESTRICT y, |
2705 | | const double * rowScale, |
2706 | | const double * columnScale) const |
2707 | | { |
2708 | | printf("createMatrix at file %s line %d\n",__FILE__,__LINE__); |
2709 | | createMatrix()->times(scalar,x,y); |
2710 | | } |
2711 | | void |
2712 | | ClpPoolMatrix::transposeTimes( double scalar, |
2713 | | const double * COIN_RESTRICT x, double * COIN_RESTRICT y, |
2714 | | const double * rowScale, |
2715 | | const double * columnScale, |
2716 | | double * spare) const |
| 2666 | void ClpPoolMatrix::times(double scalar, |
| 2667 | const double *COIN_RESTRICT x, double *COIN_RESTRICT y, |
| 2668 | const double *rowScale, |
| 2669 | const double *columnScale) const |
| 2670 | { |
| 2671 | printf("createMatrix at file %s line %d\n", __FILE__, __LINE__); |
| 2672 | createMatrix()->times(scalar, x, y); |
| 2673 | } |
| 2674 | void ClpPoolMatrix::transposeTimes(double scalar, |
| 2675 | const double *COIN_RESTRICT x, double *COIN_RESTRICT y, |
| 2676 | const double *rowScale, |
| 2677 | const double *columnScale, |
| 2678 | double *spare) const |
2904 | | void |
2905 | | ClpPoolMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar, |
2906 | | const CoinIndexedVector * rowArray, |
2907 | | CoinIndexedVector * y, |
2908 | | CoinIndexedVector * columnArray) const |
2909 | | { |
2910 | | printf("createMatrix at file %s line %d\n",__FILE__,__LINE__); |
2911 | | createMatrix()->transposeTimesByRow(model,scalar,rowArray,y,columnArray); |
| 2865 | void ClpPoolMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar, |
| 2866 | const CoinIndexedVector *rowArray, |
| 2867 | CoinIndexedVector *y, |
| 2868 | CoinIndexedVector *columnArray) const |
| 2869 | { |
| 2870 | printf("createMatrix at file %s line %d\n", __FILE__, __LINE__); |
| 2871 | createMatrix()->transposeTimesByRow(model, scalar, rowArray, y, columnArray); |
3291 | | void |
3292 | | ClpPoolMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction, |
3293 | | int & bestSequence, int & numberWanted) |
3294 | | { |
3295 | | printf("createMatrix at file %s line %d\n",__FILE__,__LINE__); |
3296 | | createMatrix()->partialPricing(model,startFraction, |
3297 | | endFraction,bestSequence, |
3298 | | numberWanted); |
| 3240 | void ClpPoolMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction, |
| 3241 | int &bestSequence, int &numberWanted) |
| 3242 | { |
| 3243 | printf("createMatrix at file %s line %d\n", __FILE__, __LINE__); |
| 3244 | createMatrix()->partialPricing(model, startFraction, |
| 3245 | endFraction, bestSequence, |
| 3246 | numberWanted); |
3312 | | int |
3313 | | ClpPoolMatrix::transposeTimes2(const ClpSimplex * model, |
3314 | | const CoinIndexedVector * pi1, CoinIndexedVector * dj1, |
3315 | | const CoinIndexedVector * pi2, |
3316 | | CoinIndexedVector * spare, |
3317 | | double * infeas, double * reducedCost, |
3318 | | double referenceIn, double devex, |
3319 | | // Array for exact devex to say what is in reference framework |
3320 | | unsigned int * COIN_RESTRICT reference, |
3321 | | double * COIN_RESTRICT weights, double scaleFactor) |
| 3259 | int ClpPoolMatrix::transposeTimes2(const ClpSimplex *model, |
| 3260 | const CoinIndexedVector *pi1, CoinIndexedVector *dj1, |
| 3261 | const CoinIndexedVector *pi2, |
| 3262 | CoinIndexedVector *spare, |
| 3263 | double *infeas, double *reducedCost, |
| 3264 | double referenceIn, double devex, |
| 3265 | // Array for exact devex to say what is in reference framework |
| 3266 | unsigned int *COIN_RESTRICT reference, |
| 3267 | double *COIN_RESTRICT weights, double scaleFactor) |
3382 | | // and do other array |
3383 | | double modification = 0.0; |
3384 | | for (j = start; j < end; j++) { |
3385 | | int iRow = stuff_[j].row_; |
3386 | | modification += piWeight[iRow]*elements_[stuff_[j].pool_]; |
3387 | | } |
3388 | | double thisWeight = weights[iColumn]; |
3389 | | double pivot = value * scaleFactor; |
3390 | | double pivotSquared = pivot * pivot; |
3391 | | thisWeight += pivotSquared * devex + pivot * modification; |
3392 | | if (thisWeight < DEVEX_TRY_NORM) { |
3393 | | if (referenceIn < 0.0) { |
3394 | | // steepest |
3395 | | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
3396 | | } else { |
3397 | | // exact |
3398 | | thisWeight = referenceIn * pivotSquared; |
3399 | | if (reference(iColumn)) |
3400 | | thisWeight += 1.0; |
3401 | | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
3402 | | } |
3403 | | } |
3404 | | weights[iColumn] = thisWeight; |
3405 | | if (!killDjs) { |
3406 | | value = reducedCost[iColumn]-value; |
3407 | | reducedCost[iColumn] = value; |
3408 | | // simplify status |
3409 | | ClpSimplex::Status status = model->getStatus(iColumn); |
3410 | | |
3411 | | switch(status) { |
3412 | | |
3413 | | case ClpSimplex::basic: |
3414 | | case ClpSimplex::isFixed: |
3415 | | break; |
3416 | | case ClpSimplex::isFree: |
3417 | | case ClpSimplex::superBasic: |
3418 | | if (fabs(value) > FREE_ACCEPT * dualTolerance) { |
3419 | | // we are going to bias towards free (but only if reasonable) |
3420 | | value *= FREE_BIAS; |
3421 | | value *= value; |
3422 | | // store square in list |
3423 | | if (infeas[iColumn]) { |
3424 | | infeas[iColumn] = value; // already there |
3425 | | } else { |
3426 | | array[numberNonZero]=value; |
3427 | | index[numberNonZero++]=iColumn; |
3428 | | } |
3429 | | } else { |
3430 | | array[numberNonZero]=0.0; |
3431 | | index[numberNonZero++]=iColumn; |
3432 | | } |
3433 | | break; |
3434 | | case ClpSimplex::atUpperBound: |
3435 | | if (value > dualTolerance) { |
3436 | | value *= value; |
3437 | | // store square in list |
3438 | | if (infeas[iColumn]) { |
3439 | | infeas[iColumn] = value; // already there |
3440 | | } else { |
3441 | | array[numberNonZero]=value; |
3442 | | index[numberNonZero++]=iColumn; |
3443 | | } |
3444 | | } else { |
3445 | | array[numberNonZero]=0.0; |
3446 | | index[numberNonZero++]=iColumn; |
3447 | | } |
3448 | | break; |
3449 | | case ClpSimplex::atLowerBound: |
3450 | | if (value < -dualTolerance) { |
3451 | | value *= value; |
3452 | | // store square in list |
3453 | | if (infeas[iColumn]) { |
3454 | | infeas[iColumn] = value; // already there |
3455 | | } else { |
3456 | | array[numberNonZero]=value; |
3457 | | index[numberNonZero++]=iColumn; |
3458 | | } |
3459 | | } else { |
3460 | | array[numberNonZero]=0.0; |
3461 | | index[numberNonZero++]=iColumn; |
3462 | | } |
3463 | | } |
3464 | | } |
| 3329 | // and do other array |
| 3330 | double modification = 0.0; |
| 3331 | for (j = start; j < end; j++) { |
| 3332 | int iRow = stuff_[j].row_; |
| 3333 | modification += piWeight[iRow] * elements_[stuff_[j].pool_]; |
| 3334 | } |
| 3335 | double thisWeight = weights[iColumn]; |
| 3336 | double pivot = value * scaleFactor; |
| 3337 | double pivotSquared = pivot * pivot; |
| 3338 | thisWeight += pivotSquared * devex + pivot * modification; |
| 3339 | if (thisWeight < DEVEX_TRY_NORM) { |
| 3340 | if (referenceIn < 0.0) { |
| 3341 | // steepest |
| 3342 | thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared); |
| 3343 | } else { |
| 3344 | // exact |
| 3345 | thisWeight = referenceIn * pivotSquared; |
| 3346 | if (reference(iColumn)) |
| 3347 | thisWeight += 1.0; |
| 3348 | thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM); |
| 3349 | } |
| 3350 | } |
| 3351 | weights[iColumn] = thisWeight; |
| 3352 | if (!killDjs) { |
| 3353 | value = reducedCost[iColumn] - value; |
| 3354 | reducedCost[iColumn] = value; |
| 3355 | // simplify status |
| 3356 | ClpSimplex::Status status = model->getStatus(iColumn); |
| 3357 | |
| 3358 | switch (status) { |
| 3359 | |
| 3360 | case ClpSimplex::basic: |
| 3361 | case ClpSimplex::isFixed: |
| 3362 | break; |
| 3363 | case ClpSimplex::isFree: |
| 3364 | case ClpSimplex::superBasic: |
| 3365 | if (fabs(value) > FREE_ACCEPT * dualTolerance) { |
| 3366 | // we are going to bias towards free (but only if reasonable) |
| 3367 | value *= FREE_BIAS; |
| 3368 | value *= value; |
| 3369 | // store square in list |
| 3370 | if (infeas[iColumn]) { |
| 3371 | infeas[iColumn] = value; // already there |
| 3372 | } else { |
| 3373 | array[numberNonZero] = value; |
| 3374 | index[numberNonZero++] = iColumn; |
| 3375 | } |
| 3376 | } else { |
| 3377 | array[numberNonZero] = 0.0; |
| 3378 | index[numberNonZero++] = iColumn; |
| 3379 | } |
| 3380 | break; |
| 3381 | case ClpSimplex::atUpperBound: |
| 3382 | if (value > dualTolerance) { |
| 3383 | value *= value; |
| 3384 | // store square in list |
| 3385 | if (infeas[iColumn]) { |
| 3386 | infeas[iColumn] = value; // already there |
| 3387 | } else { |
| 3388 | array[numberNonZero] = value; |
| 3389 | index[numberNonZero++] = iColumn; |
| 3390 | } |
| 3391 | } else { |
| 3392 | array[numberNonZero] = 0.0; |
| 3393 | index[numberNonZero++] = iColumn; |
| 3394 | } |
| 3395 | break; |
| 3396 | case ClpSimplex::atLowerBound: |
| 3397 | if (value < -dualTolerance) { |
| 3398 | value *= value; |
| 3399 | // store square in list |
| 3400 | if (infeas[iColumn]) { |
| 3401 | |