1 | // Copyright (C) 2003, International Business Machines |
---|
2 | // Corporation and others. All Rights Reserved. |
---|
3 | |
---|
4 | |
---|
5 | #include <cstdio> |
---|
6 | |
---|
7 | #include "CoinPragma.hpp" |
---|
8 | #include "CoinIndexedVector.hpp" |
---|
9 | #include "CoinHelperFunctions.hpp" |
---|
10 | |
---|
11 | #include "ClpSimplex.hpp" |
---|
12 | #include "ClpFactorization.hpp" |
---|
13 | // at end to get min/max! |
---|
14 | #include "ClpPlusMinusOneMatrix.hpp" |
---|
15 | #include "ClpMessage.hpp" |
---|
16 | |
---|
17 | //############################################################################# |
---|
18 | // Constructors / Destructor / Assignment |
---|
19 | //############################################################################# |
---|
20 | |
---|
21 | //------------------------------------------------------------------- |
---|
22 | // Default Constructor |
---|
23 | //------------------------------------------------------------------- |
---|
24 | ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix () |
---|
25 | : ClpMatrixBase() |
---|
26 | { |
---|
27 | setType(12); |
---|
28 | elements_ = NULL; |
---|
29 | startPositive_ = NULL; |
---|
30 | startNegative_ = NULL; |
---|
31 | lengths_=NULL; |
---|
32 | indices_=NULL; |
---|
33 | numberRows_=0; |
---|
34 | numberColumns_=0; |
---|
35 | columnOrdered_=true; |
---|
36 | } |
---|
37 | |
---|
38 | //------------------------------------------------------------------- |
---|
39 | // Copy constructor |
---|
40 | //------------------------------------------------------------------- |
---|
41 | ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const ClpPlusMinusOneMatrix & rhs) |
---|
42 | : ClpMatrixBase(rhs) |
---|
43 | { |
---|
44 | elements_ = NULL; |
---|
45 | startPositive_ = NULL; |
---|
46 | startNegative_ = NULL; |
---|
47 | lengths_=NULL; |
---|
48 | indices_=NULL; |
---|
49 | numberRows_=rhs.numberRows_; |
---|
50 | numberColumns_=rhs.numberColumns_; |
---|
51 | columnOrdered_=rhs.columnOrdered_; |
---|
52 | if (numberColumns_) { |
---|
53 | indices_ = new int [ 2*numberColumns_]; |
---|
54 | memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int)); |
---|
55 | startPositive_ = new int [ numberColumns_+1]; |
---|
56 | memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int)); |
---|
57 | startNegative_ = new int [ numberColumns_]; |
---|
58 | memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(int)); |
---|
59 | } |
---|
60 | } |
---|
61 | |
---|
62 | ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const CoinPackedMatrix & rhs) |
---|
63 | : ClpMatrixBase() |
---|
64 | { |
---|
65 | setType(12); |
---|
66 | elements_ = NULL; |
---|
67 | startPositive_ = NULL; |
---|
68 | startNegative_ = NULL; |
---|
69 | lengths_=NULL; |
---|
70 | indices_=NULL; |
---|
71 | int iColumn; |
---|
72 | assert (rhs.isColOrdered()); |
---|
73 | // get matrix data pointers |
---|
74 | const int * row = rhs.getIndices(); |
---|
75 | const CoinBigIndex * columnStart = rhs.getVectorStarts(); |
---|
76 | const int * columnLength = rhs.getVectorLengths(); |
---|
77 | const double * elementByColumn = rhs.getElements(); |
---|
78 | numberColumns_ = rhs.getNumCols(); |
---|
79 | bool goodPlusMinusOne=true; |
---|
80 | numberRows_=-1; |
---|
81 | indices_ = new int[rhs.getNumElements()]; |
---|
82 | startPositive_ = new int [numberColumns_+1]; |
---|
83 | startNegative_ = new int [numberColumns_]; |
---|
84 | int * temp = new int [rhs.getNumRows()]; |
---|
85 | CoinBigIndex j=0; |
---|
86 | for (iColumn=0;iColumn<numberColumns_;iColumn++) { |
---|
87 | CoinBigIndex k; |
---|
88 | int iNeg=0; |
---|
89 | startPositive_[iColumn]=j; |
---|
90 | for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn]; |
---|
91 | k++) { |
---|
92 | int iRow; |
---|
93 | if (fabs(elementByColumn[k]-1.0)<1.0e-10) { |
---|
94 | iRow = row[k]; |
---|
95 | numberRows_ = max(numberRows_,iRow); |
---|
96 | indices_[j++]=iRow; |
---|
97 | } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) { |
---|
98 | iRow = row[k]; |
---|
99 | numberRows_ = max(numberRows_,iRow); |
---|
100 | temp[iNeg++]=iRow; |
---|
101 | } else { |
---|
102 | goodPlusMinusOne = false; // not a network |
---|
103 | } |
---|
104 | } |
---|
105 | if (goodPlusMinusOne) { |
---|
106 | // move negative |
---|
107 | startNegative_[iColumn]=j; |
---|
108 | for (k=0;k<iNeg;k++) { |
---|
109 | indices_[j++] = temp[k]; |
---|
110 | } |
---|
111 | } else { |
---|
112 | break; |
---|
113 | } |
---|
114 | } |
---|
115 | startPositive_[numberColumns_]=j; |
---|
116 | delete [] temp; |
---|
117 | if (!goodPlusMinusOne) { |
---|
118 | delete [] indices_; |
---|
119 | // put in message |
---|
120 | printf("Not all +-1 - can test if indices_ null\n"); |
---|
121 | indices_=NULL; |
---|
122 | numberRows_=0; |
---|
123 | numberColumns_=0; |
---|
124 | delete [] startPositive_; |
---|
125 | delete [] startNegative_; |
---|
126 | startPositive_ = NULL; |
---|
127 | startNegative_ = NULL; |
---|
128 | } else { |
---|
129 | numberRows_ ++; // correct |
---|
130 | columnOrdered_ = true; |
---|
131 | } |
---|
132 | } |
---|
133 | |
---|
134 | //------------------------------------------------------------------- |
---|
135 | // Destructor |
---|
136 | //------------------------------------------------------------------- |
---|
137 | ClpPlusMinusOneMatrix::~ClpPlusMinusOneMatrix () |
---|
138 | { |
---|
139 | delete [] elements_; |
---|
140 | delete [] startPositive_; |
---|
141 | delete [] startNegative_; |
---|
142 | delete [] lengths_; |
---|
143 | delete [] indices_; |
---|
144 | } |
---|
145 | |
---|
146 | //---------------------------------------------------------------- |
---|
147 | // Assignment operator |
---|
148 | //------------------------------------------------------------------- |
---|
149 | ClpPlusMinusOneMatrix & |
---|
150 | ClpPlusMinusOneMatrix::operator=(const ClpPlusMinusOneMatrix& rhs) |
---|
151 | { |
---|
152 | if (this != &rhs) { |
---|
153 | ClpMatrixBase::operator=(rhs); |
---|
154 | delete [] elements_; |
---|
155 | delete [] startPositive_; |
---|
156 | delete [] startNegative_; |
---|
157 | delete [] lengths_; |
---|
158 | delete [] indices_; |
---|
159 | elements_ = NULL; |
---|
160 | startPositive_ = NULL; |
---|
161 | lengths_=NULL; |
---|
162 | indices_=NULL; |
---|
163 | numberRows_=rhs.numberRows_; |
---|
164 | numberColumns_=rhs.numberColumns_; |
---|
165 | columnOrdered_=rhs.columnOrdered_; |
---|
166 | if (numberColumns_) { |
---|
167 | indices_ = new int [ 2*numberColumns_]; |
---|
168 | memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int)); |
---|
169 | startPositive_ = new int [ numberColumns_+1]; |
---|
170 | memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int)); |
---|
171 | startNegative_ = new int [ numberColumns_]; |
---|
172 | memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(int)); |
---|
173 | } |
---|
174 | } |
---|
175 | return *this; |
---|
176 | } |
---|
177 | //------------------------------------------------------------------- |
---|
178 | // Clone |
---|
179 | //------------------------------------------------------------------- |
---|
180 | ClpMatrixBase * ClpPlusMinusOneMatrix::clone() const |
---|
181 | { |
---|
182 | return new ClpPlusMinusOneMatrix(*this); |
---|
183 | } |
---|
184 | |
---|
185 | /* Returns a new matrix in reverse order without gaps */ |
---|
186 | ClpMatrixBase * |
---|
187 | ClpPlusMinusOneMatrix::reverseOrderedCopy() const |
---|
188 | { |
---|
189 | int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_; |
---|
190 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
---|
191 | // count number in each row/column |
---|
192 | int * tempP = new int [numberMinor]; |
---|
193 | int * tempN = new int [numberMinor]; |
---|
194 | memset(tempP,0,numberMinor*sizeof(int)); |
---|
195 | memset(tempN,0,numberMinor*sizeof(int)); |
---|
196 | CoinBigIndex j=0; |
---|
197 | int i; |
---|
198 | for (i=0;i<numberMajor;i++) { |
---|
199 | for (;j<startNegative_[i];j++) { |
---|
200 | int iRow = indices_[j]; |
---|
201 | tempP[iRow]++; |
---|
202 | } |
---|
203 | for (;j<startPositive_[i+1];j++) { |
---|
204 | int iRow = indices_[j]; |
---|
205 | tempN[iRow]++; |
---|
206 | } |
---|
207 | } |
---|
208 | int * newIndices = new int [startPositive_[numberMajor]]; |
---|
209 | int * newP = new int [numberMinor+1]; |
---|
210 | int * newN = new int[numberMinor]; |
---|
211 | int iRow; |
---|
212 | j=0; |
---|
213 | // do starts |
---|
214 | for (iRow=0;iRow<numberMinor;iRow++) { |
---|
215 | newP[iRow]=j; |
---|
216 | j += tempP[iRow]; |
---|
217 | tempP[iRow]=newP[iRow]; |
---|
218 | newN[iRow] = j; |
---|
219 | j += tempN[iRow]; |
---|
220 | tempN[iRow]=newN[iRow]; |
---|
221 | } |
---|
222 | newP[numberMinor]=j; |
---|
223 | j=0; |
---|
224 | for (i=0;i<numberMajor;i++) { |
---|
225 | for (;j<startNegative_[i];j++) { |
---|
226 | int iRow = indices_[j]; |
---|
227 | int put = tempP[iRow]; |
---|
228 | newIndices[put++] = i; |
---|
229 | tempP[iRow] = put; |
---|
230 | } |
---|
231 | for (;j<startPositive_[i+1];j++) { |
---|
232 | int iRow = indices_[j]; |
---|
233 | int put = tempN[iRow]; |
---|
234 | newIndices[put++] = i; |
---|
235 | tempN[iRow] = put; |
---|
236 | } |
---|
237 | } |
---|
238 | delete [] tempP; |
---|
239 | delete [] tempN; |
---|
240 | ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix(); |
---|
241 | newCopy->passInCopy(numberMinor, numberMajor, |
---|
242 | !columnOrdered_, newIndices, newP, newN); |
---|
243 | return newCopy; |
---|
244 | } |
---|
245 | //unscaled versions |
---|
246 | void |
---|
247 | ClpPlusMinusOneMatrix::times(double scalar, |
---|
248 | const double * x, double * y) const |
---|
249 | { |
---|
250 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
---|
251 | int i; |
---|
252 | CoinBigIndex j; |
---|
253 | assert (columnOrdered_); |
---|
254 | for (i=0;i<numberMajor;i++) { |
---|
255 | double value = scalar*x[i]; |
---|
256 | if (value) { |
---|
257 | for (j=startPositive_[i];j<startNegative_[i];j++) { |
---|
258 | int iRow = indices_[j]; |
---|
259 | y[iRow] += value; |
---|
260 | } |
---|
261 | for (;j<startPositive_[i+1];j++) { |
---|
262 | int iRow = indices_[j]; |
---|
263 | y[iRow] -= value; |
---|
264 | } |
---|
265 | } |
---|
266 | } |
---|
267 | } |
---|
268 | void |
---|
269 | ClpPlusMinusOneMatrix::transposeTimes(double scalar, |
---|
270 | const double * x, double * y) const |
---|
271 | { |
---|
272 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
---|
273 | int i; |
---|
274 | CoinBigIndex j=0; |
---|
275 | assert (columnOrdered_); |
---|
276 | for (i=0;i<numberMajor;i++) { |
---|
277 | double value = 0.0; |
---|
278 | for (;j<startNegative_[i];j++) { |
---|
279 | int iRow = indices_[j]; |
---|
280 | value += x[iRow]; |
---|
281 | } |
---|
282 | for (;j<startPositive_[i+1];j++) { |
---|
283 | int iRow = indices_[j]; |
---|
284 | value -= x[iRow]; |
---|
285 | } |
---|
286 | y[i] += scalar*value; |
---|
287 | } |
---|
288 | } |
---|
289 | void |
---|
290 | ClpPlusMinusOneMatrix::times(double scalar, |
---|
291 | const double * x, double * y, |
---|
292 | const double * rowScale, |
---|
293 | const double * columnScale) const |
---|
294 | { |
---|
295 | // we know it is not scaled |
---|
296 | times(scalar, x, y); |
---|
297 | } |
---|
298 | void |
---|
299 | ClpPlusMinusOneMatrix::transposeTimes( double scalar, |
---|
300 | const double * x, double * y, |
---|
301 | const double * rowScale, |
---|
302 | const double * columnScale) const |
---|
303 | { |
---|
304 | // we know it is not scaled |
---|
305 | transposeTimes(scalar, x, y); |
---|
306 | } |
---|
307 | /* Return <code>x * A + y</code> in <code>z</code>. |
---|
308 | Squashes small elements and knows about ClpSimplex */ |
---|
309 | void |
---|
310 | ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex * model, double scalar, |
---|
311 | const CoinIndexedVector * rowArray, |
---|
312 | CoinIndexedVector * y, |
---|
313 | CoinIndexedVector * columnArray) const |
---|
314 | { |
---|
315 | // we know it is not scaled |
---|
316 | columnArray->clear(); |
---|
317 | double * pi = rowArray->denseVector(); |
---|
318 | int numberNonZero=0; |
---|
319 | int * index = columnArray->getIndices(); |
---|
320 | double * array = columnArray->denseVector(); |
---|
321 | int numberInRowArray = rowArray->getNumElements(); |
---|
322 | // maybe I need one in OsiSimplex |
---|
323 | double zeroTolerance = model->factorization()->zeroTolerance(); |
---|
324 | int numberRows = model->numberRows(); |
---|
325 | ClpPlusMinusOneMatrix* rowCopy = |
---|
326 | dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy()); |
---|
327 | if (numberInRowArray>0.3*numberRows||!rowCopy) { |
---|
328 | // do by column |
---|
329 | int iColumn; |
---|
330 | double * markVector = y->denseVector(); // probably empty |
---|
331 | CoinBigIndex j=0; |
---|
332 | assert (columnOrdered_); |
---|
333 | for (iColumn=0;iColumn<numberColumns_;iColumn++) { |
---|
334 | double value2 = 0.0; |
---|
335 | for (;j<startNegative_[iColumn];j++) { |
---|
336 | int iRow = indices_[j]; |
---|
337 | value2 += pi[iRow]; |
---|
338 | } |
---|
339 | for (;j<startPositive_[iColumn+1];j++) { |
---|
340 | int iRow = indices_[j]; |
---|
341 | value2 -= pi[iRow]; |
---|
342 | } |
---|
343 | double value = markVector[iColumn]; |
---|
344 | markVector[iColumn]=0.0; |
---|
345 | value += scalar*value2; |
---|
346 | if (fabs(value)>zeroTolerance) { |
---|
347 | index[numberNonZero++]=iColumn; |
---|
348 | array[iColumn]=value; |
---|
349 | } |
---|
350 | } |
---|
351 | columnArray->setNumElements(numberNonZero); |
---|
352 | y->setNumElements(0); |
---|
353 | } else { |
---|
354 | // do by row |
---|
355 | rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray); |
---|
356 | } |
---|
357 | } |
---|
358 | /* Return <code>x * A + y</code> in <code>z</code>. |
---|
359 | Squashes small elements and knows about ClpSimplex */ |
---|
360 | void |
---|
361 | ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar, |
---|
362 | const CoinIndexedVector * rowArray, |
---|
363 | CoinIndexedVector * y, |
---|
364 | CoinIndexedVector * columnArray) const |
---|
365 | { |
---|
366 | columnArray->clear(); |
---|
367 | double * pi = rowArray->denseVector(); |
---|
368 | int numberNonZero=0; |
---|
369 | int * index = columnArray->getIndices(); |
---|
370 | double * array = columnArray->denseVector(); |
---|
371 | int numberInRowArray = rowArray->getNumElements(); |
---|
372 | // maybe I need one in OsiSimplex |
---|
373 | double zeroTolerance = model->factorization()->zeroTolerance(); |
---|
374 | const int * column = indices_; |
---|
375 | const CoinBigIndex * startPositive = startPositive_; |
---|
376 | const CoinBigIndex * startNegative = startNegative_; |
---|
377 | const int * whichRow = rowArray->getIndices(); |
---|
378 | if (numberInRowArray>2||y->getNumElements()) { |
---|
379 | // do by rows |
---|
380 | int iRow; |
---|
381 | double * markVector = y->denseVector(); // probably empty .. but |
---|
382 | int * mark = y->getIndices(); |
---|
383 | int numberOriginal=y->getNumElements(); |
---|
384 | int i; |
---|
385 | for (i=0;i<numberOriginal;i++) { |
---|
386 | int iColumn = mark[i]; |
---|
387 | index[i]=iColumn; |
---|
388 | array[iColumn]=markVector[iColumn]; |
---|
389 | markVector[iColumn]=0.0; |
---|
390 | } |
---|
391 | numberNonZero=numberOriginal; |
---|
392 | // and set up mark as char array |
---|
393 | char * marked = (char *) markVector; |
---|
394 | for (i=0;i<numberOriginal;i++) { |
---|
395 | int iColumn = index[i]; |
---|
396 | marked[iColumn]=0; |
---|
397 | } |
---|
398 | for (i=0;i<numberInRowArray;i++) { |
---|
399 | iRow = whichRow[i]; |
---|
400 | double value = pi[iRow]*scalar; |
---|
401 | CoinBigIndex j; |
---|
402 | for (j=startPositive[iRow];j<startNegative[iRow];j++) { |
---|
403 | int iColumn = column[j]; |
---|
404 | if (!marked[iColumn]) { |
---|
405 | marked[iColumn]=1; |
---|
406 | index[numberNonZero++]=iColumn; |
---|
407 | } |
---|
408 | array[iColumn] += value; |
---|
409 | } |
---|
410 | for (j=startNegative[iRow];j<startPositive[iRow+1];j++) { |
---|
411 | int iColumn = column[j]; |
---|
412 | if (!marked[iColumn]) { |
---|
413 | marked[iColumn]=1; |
---|
414 | index[numberNonZero++]=iColumn; |
---|
415 | } |
---|
416 | array[iColumn] -= value; |
---|
417 | } |
---|
418 | } |
---|
419 | // get rid of tiny values and zero out marked |
---|
420 | numberOriginal=numberNonZero; |
---|
421 | numberNonZero=0; |
---|
422 | for (i=0;i<numberOriginal;i++) { |
---|
423 | int iColumn = index[i]; |
---|
424 | marked[iColumn]=0; |
---|
425 | if (fabs(array[iColumn])>zeroTolerance) { |
---|
426 | index[numberNonZero++]=iColumn; |
---|
427 | } else { |
---|
428 | array[iColumn]=0.0; |
---|
429 | } |
---|
430 | } |
---|
431 | } else if (numberInRowArray==2) { |
---|
432 | // do by rows when two rows (do longer first) |
---|
433 | int iRow0 = whichRow[0]; |
---|
434 | int iRow1 = whichRow[1]; |
---|
435 | if (startPositive[iRow0+1]-startPositive[iRow0]< |
---|
436 | startPositive[iRow1+1]-startPositive[iRow1]) { |
---|
437 | int temp = iRow0; |
---|
438 | iRow0 = iRow1; |
---|
439 | iRow1 = temp; |
---|
440 | } |
---|
441 | int numberOriginal; |
---|
442 | int i; |
---|
443 | numberNonZero=0; |
---|
444 | double value; |
---|
445 | value = pi[iRow0]*scalar; |
---|
446 | CoinBigIndex j; |
---|
447 | for (j=startPositive[iRow0];j<startNegative[iRow0];j++) { |
---|
448 | int iColumn = column[j]; |
---|
449 | index[numberNonZero++]=iColumn; |
---|
450 | array[iColumn] = value; |
---|
451 | } |
---|
452 | for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) { |
---|
453 | int iColumn = column[j]; |
---|
454 | index[numberNonZero++]=iColumn; |
---|
455 | array[iColumn] = -value; |
---|
456 | } |
---|
457 | value = pi[iRow1]*scalar; |
---|
458 | for (j=startPositive[iRow1];j<startNegative[iRow1];j++) { |
---|
459 | int iColumn = column[j]; |
---|
460 | double value2= array[iColumn]; |
---|
461 | if (value2) { |
---|
462 | value2 += value; |
---|
463 | } else { |
---|
464 | value2 = value; |
---|
465 | index[numberNonZero++]=iColumn; |
---|
466 | } |
---|
467 | array[iColumn] = value2; |
---|
468 | } |
---|
469 | for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) { |
---|
470 | int iColumn = column[j]; |
---|
471 | double value2= array[iColumn]; |
---|
472 | if (value2) { |
---|
473 | value2 -= value; |
---|
474 | } else { |
---|
475 | value2 = -value; |
---|
476 | index[numberNonZero++]=iColumn; |
---|
477 | } |
---|
478 | array[iColumn] = value2; |
---|
479 | } |
---|
480 | // get rid of tiny values and zero out marked |
---|
481 | numberOriginal=numberNonZero; |
---|
482 | numberNonZero=0; |
---|
483 | for (i=0;i<numberOriginal;i++) { |
---|
484 | int iColumn = index[i]; |
---|
485 | if (fabs(array[iColumn])>zeroTolerance) { |
---|
486 | index[numberNonZero++]=iColumn; |
---|
487 | } else { |
---|
488 | array[iColumn]=0.0; |
---|
489 | } |
---|
490 | } |
---|
491 | } else if (numberInRowArray==1) { |
---|
492 | // Just one row |
---|
493 | int iRow=rowArray->getIndices()[0]; |
---|
494 | numberNonZero=0; |
---|
495 | double value; |
---|
496 | iRow = whichRow[0]; |
---|
497 | value = pi[iRow]*scalar; |
---|
498 | if (fabs(value)>zeroTolerance) { |
---|
499 | CoinBigIndex j; |
---|
500 | for (j=startPositive[iRow];j<startNegative[iRow];j++) { |
---|
501 | int iColumn = column[j]; |
---|
502 | index[numberNonZero++]=iColumn; |
---|
503 | array[iColumn] = value; |
---|
504 | } |
---|
505 | for (j=startNegative[iRow];j<startPositive[iRow+1];j++) { |
---|
506 | int iColumn = column[j]; |
---|
507 | index[numberNonZero++]=iColumn; |
---|
508 | array[iColumn] = -value; |
---|
509 | } |
---|
510 | } |
---|
511 | } |
---|
512 | columnArray->setNumElements(numberNonZero); |
---|
513 | y->setNumElements(0); |
---|
514 | } |
---|
515 | /* Return <code>x *A in <code>z</code> but |
---|
516 | just for indices in y. |
---|
517 | Squashes small elements and knows about ClpSimplex */ |
---|
518 | void |
---|
519 | ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex * model, |
---|
520 | const CoinIndexedVector * rowArray, |
---|
521 | const CoinIndexedVector * y, |
---|
522 | CoinIndexedVector * columnArray) const |
---|
523 | { |
---|
524 | columnArray->clear(); |
---|
525 | double * pi = rowArray->denseVector(); |
---|
526 | int numberNonZero=0; |
---|
527 | int * index = columnArray->getIndices(); |
---|
528 | double * array = columnArray->denseVector(); |
---|
529 | // maybe I need one in OsiSimplex |
---|
530 | double zeroTolerance = model->factorization()->zeroTolerance(); |
---|
531 | int jColumn; |
---|
532 | int numberToDo = y->getNumElements(); |
---|
533 | const int * which = y->getIndices(); |
---|
534 | for (jColumn=0;jColumn<numberToDo;jColumn++) { |
---|
535 | int iColumn = which[jColumn]; |
---|
536 | double value = 0.0; |
---|
537 | CoinBigIndex j=startPositive_[iColumn]; |
---|
538 | for (;j<startNegative_[iColumn];j++) { |
---|
539 | int iRow = indices_[j]; |
---|
540 | value += pi[iRow]; |
---|
541 | } |
---|
542 | for (;j<startPositive_[iColumn+1];j++) { |
---|
543 | int iRow = indices_[j]; |
---|
544 | value -= pi[iRow]; |
---|
545 | } |
---|
546 | if (fabs(value)>zeroTolerance) { |
---|
547 | index[numberNonZero++]=iColumn; |
---|
548 | array[iColumn]=value; |
---|
549 | } |
---|
550 | } |
---|
551 | } |
---|
552 | /* Returns number of elements in basis |
---|
553 | column is basic if entry >=0 */ |
---|
554 | CoinBigIndex |
---|
555 | ClpPlusMinusOneMatrix::numberInBasis(const int * columnIsBasic) const |
---|
556 | { |
---|
557 | int i; |
---|
558 | CoinBigIndex numberElements=0; |
---|
559 | assert (columnOrdered_); |
---|
560 | for (i=0;i<numberColumns_;i++) { |
---|
561 | if (columnIsBasic[i]>=0) |
---|
562 | numberElements += startPositive_[i+1]-startPositive_[i]; |
---|
563 | } |
---|
564 | return numberElements; |
---|
565 | } |
---|
566 | // Fills in basis (Returns number of elements and updates numberBasic) |
---|
567 | CoinBigIndex |
---|
568 | ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model, |
---|
569 | const int * columnIsBasic, int & numberBasic, |
---|
570 | int * indexRowU, int * indexColumnU, |
---|
571 | double * elementU) const |
---|
572 | { |
---|
573 | #ifdef CLPDEBUG |
---|
574 | const double * rowScale = model->rowScale(); |
---|
575 | assert (!rowScale); |
---|
576 | #endif |
---|
577 | int i; |
---|
578 | CoinBigIndex numberElements=0; |
---|
579 | assert (columnOrdered_); |
---|
580 | for (i=0;i<numberColumns_;i++) { |
---|
581 | if (columnIsBasic[i]>=0) { |
---|
582 | CoinBigIndex j=startPositive_[i]; |
---|
583 | for (;j<startNegative_[i];j++) { |
---|
584 | int iRow = indices_[j]; |
---|
585 | indexRowU[numberElements]=iRow; |
---|
586 | indexColumnU[numberElements]=numberBasic; |
---|
587 | elementU[numberElements++]=1.0; |
---|
588 | } |
---|
589 | for (;j<startPositive_[i+1];j++) { |
---|
590 | int iRow = indices_[j]; |
---|
591 | indexRowU[numberElements]=iRow; |
---|
592 | indexColumnU[numberElements]=numberBasic; |
---|
593 | elementU[numberElements++]=-1.0; |
---|
594 | } |
---|
595 | numberBasic++; |
---|
596 | } |
---|
597 | } |
---|
598 | return numberElements; |
---|
599 | } |
---|
600 | /* Unpacks a column into an CoinIndexedvector |
---|
601 | Note that model is NOT const. Bounds and objective could |
---|
602 | be modified if doing column generation */ |
---|
603 | void |
---|
604 | ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model, |
---|
605 | CoinIndexedVector * rowArray, |
---|
606 | int iColumn) const |
---|
607 | { |
---|
608 | CoinBigIndex j=startPositive_[iColumn]; |
---|
609 | for (;j<startNegative_[iColumn];j++) { |
---|
610 | int iRow = indices_[j]; |
---|
611 | rowArray->add(iRow,1.0); |
---|
612 | } |
---|
613 | for (;j<startPositive_[iColumn+1];j++) { |
---|
614 | int iRow = indices_[j]; |
---|
615 | rowArray->add(iRow,-1.0); |
---|
616 | } |
---|
617 | } |
---|
618 | /* Adds multiple of a column into an CoinIndexedvector |
---|
619 | You can use quickAdd to add to vector */ |
---|
620 | void |
---|
621 | ClpPlusMinusOneMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray, |
---|
622 | int iColumn, double multiplier) const |
---|
623 | { |
---|
624 | CoinBigIndex j=startPositive_[iColumn]; |
---|
625 | for (;j<startNegative_[iColumn];j++) { |
---|
626 | int iRow = indices_[j]; |
---|
627 | rowArray->quickAdd(iRow,multiplier); |
---|
628 | } |
---|
629 | for (;j<startPositive_[iColumn+1];j++) { |
---|
630 | int iRow = indices_[j]; |
---|
631 | rowArray->quickAdd(iRow,-multiplier); |
---|
632 | } |
---|
633 | } |
---|
634 | |
---|
635 | // Return a complete CoinPackedMatrix |
---|
636 | CoinPackedMatrix * |
---|
637 | ClpPlusMinusOneMatrix::getPackedMatrix() const |
---|
638 | { |
---|
639 | int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_; |
---|
640 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
---|
641 | return new CoinPackedMatrix(columnOrdered_,numberMinor,numberMajor, |
---|
642 | getNumElements(), |
---|
643 | getElements(),indices_, |
---|
644 | startPositive_,getVectorLengths()); |
---|
645 | |
---|
646 | } |
---|
647 | /* A vector containing the elements in the packed matrix. Note that there |
---|
648 | might be gaps in this list, entries that do not belong to any |
---|
649 | major-dimension vector. To get the actual elements one should look at |
---|
650 | this vector together with vectorStarts and vectorLengths. */ |
---|
651 | const double * |
---|
652 | ClpPlusMinusOneMatrix::getElements() const |
---|
653 | { |
---|
654 | if (!elements_) { |
---|
655 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
---|
656 | int numberElements = startPositive_[numberMajor]; |
---|
657 | elements_ = new double [numberElements]; |
---|
658 | CoinBigIndex j=0; |
---|
659 | int i; |
---|
660 | for (i=0;i<numberMajor;i++) { |
---|
661 | for (;j<startNegative_[i];j++) { |
---|
662 | elements_[j]=1.0; |
---|
663 | } |
---|
664 | for (;j<startPositive_[i+1];j++) { |
---|
665 | elements_[j]=-1.0; |
---|
666 | } |
---|
667 | } |
---|
668 | } |
---|
669 | return elements_; |
---|
670 | } |
---|
671 | |
---|
672 | const CoinBigIndex * |
---|
673 | ClpPlusMinusOneMatrix::getVectorStarts() const |
---|
674 | { |
---|
675 | return startPositive_; |
---|
676 | } |
---|
677 | /* The lengths of the major-dimension vectors. */ |
---|
678 | const int * |
---|
679 | ClpPlusMinusOneMatrix::getVectorLengths() const |
---|
680 | { |
---|
681 | if (!lengths_) { |
---|
682 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
---|
683 | lengths_ = new int [numberMajor]; |
---|
684 | int i; |
---|
685 | for (i=0;i<numberMajor;i++) { |
---|
686 | lengths_[i]=startPositive_[i+1]-startPositive_[i]; |
---|
687 | } |
---|
688 | } |
---|
689 | return lengths_; |
---|
690 | } |
---|
691 | /* Delete the columns whose indices are listed in <code>indDel</code>. */ |
---|
692 | void |
---|
693 | ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel) |
---|
694 | { |
---|
695 | abort(); |
---|
696 | } |
---|
697 | /* Delete the rows whose indices are listed in <code>indDel</code>. */ |
---|
698 | void |
---|
699 | ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel) |
---|
700 | { |
---|
701 | abort(); |
---|
702 | } |
---|
703 | bool |
---|
704 | ClpPlusMinusOneMatrix::isColOrdered() const |
---|
705 | { |
---|
706 | return columnOrdered_; |
---|
707 | } |
---|
708 | /* Number of entries in the packed matrix. */ |
---|
709 | CoinBigIndex |
---|
710 | ClpPlusMinusOneMatrix::getNumElements() const |
---|
711 | { |
---|
712 | int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_; |
---|
713 | if (startPositive_) |
---|
714 | return startPositive_[numberMajor]; |
---|
715 | else |
---|
716 | return 0; |
---|
717 | } |
---|
718 | // pass in copy (object takes ownership) |
---|
719 | void |
---|
720 | ClpPlusMinusOneMatrix::passInCopy(int numberRows, int numberColumns, |
---|
721 | bool columnOrdered, int * indices, |
---|
722 | int * startPositive, int * startNegative) |
---|
723 | { |
---|
724 | columnOrdered_=columnOrdered; |
---|
725 | startPositive_ = startPositive; |
---|
726 | startNegative_ = startNegative; |
---|
727 | indices_ = indices; |
---|
728 | numberRows_=numberRows; |
---|
729 | numberColumns_=numberColumns; |
---|
730 | } |
---|