source: trunk/Clp/src/ClpPlusMinusOneMatrix.cpp @ 1525

Last change on this file since 1525 was 1525, checked in by mjs, 10 years ago

Formatted .cpp, .hpp, .c, .h files with "astyle -A4 -p". This matches the formatting used in the grand CBC reorganization.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 73.5 KB
Line 
1/* $Id: ClpPlusMinusOneMatrix.cpp 1525 2010-02-26 17:27:59Z mjs $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4
5
6#include <cstdio>
7
8#include "CoinPragma.hpp"
9#include "CoinIndexedVector.hpp"
10#include "CoinPackedVector.hpp"
11#include "CoinHelperFunctions.hpp"
12
13#include "ClpSimplex.hpp"
14#include "ClpFactorization.hpp"
15// at end to get min/max!
16#include "ClpPlusMinusOneMatrix.hpp"
17#include "ClpMessage.hpp"
18
19//#############################################################################
20// Constructors / Destructor / Assignment
21//#############################################################################
22
23//-------------------------------------------------------------------
24// Default Constructor
25//-------------------------------------------------------------------
26ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix ()
27     : ClpMatrixBase()
28{
29     setType(12);
30     matrix_ = NULL;
31     startPositive_ = NULL;
32     startNegative_ = NULL;
33     lengths_ = NULL;
34     indices_ = NULL;
35     numberRows_ = 0;
36     numberColumns_ = 0;
37     columnOrdered_ = true;
38}
39
40//-------------------------------------------------------------------
41// Copy constructor
42//-------------------------------------------------------------------
43ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const ClpPlusMinusOneMatrix & rhs)
44     : ClpMatrixBase(rhs)
45{
46     matrix_ = NULL;
47     startPositive_ = NULL;
48     startNegative_ = NULL;
49     lengths_ = NULL;
50     indices_ = NULL;
51     numberRows_ = rhs.numberRows_;
52     numberColumns_ = rhs.numberColumns_;
53     columnOrdered_ = rhs.columnOrdered_;
54     if (numberColumns_) {
55          CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
56          indices_ = new int [ numberElements];
57          CoinMemcpyN(rhs.indices_, numberElements, indices_);
58          startPositive_ = new CoinBigIndex [ numberColumns_+1];
59          CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_);
60          startNegative_ = new CoinBigIndex [ numberColumns_];
61          CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_);
62     }
63     int numberRows = getNumRows();
64     if (rhs.rhsOffset_ && numberRows) {
65          rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
66     } else {
67          rhsOffset_ = NULL;
68     }
69}
70// Constructor from arrays
71ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix(int numberRows, int numberColumns,
72          bool columnOrdered, const int * indices,
73          const CoinBigIndex * startPositive,
74          const CoinBigIndex * startNegative)
75     : ClpMatrixBase()
76{
77     setType(12);
78     matrix_ = NULL;
79     lengths_ = NULL;
80     numberRows_ = numberRows;
81     numberColumns_ = numberColumns;
82     columnOrdered_ = columnOrdered;
83     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
84     int numberElements = startPositive[numberMajor];
85     startPositive_ = ClpCopyOfArray(startPositive, numberMajor + 1);
86     startNegative_ = ClpCopyOfArray(startNegative, numberMajor);
87     indices_ = ClpCopyOfArray(indices, numberElements);
88     // Check valid
89     checkValid(false);
90}
91
92ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const CoinPackedMatrix & rhs)
93     : ClpMatrixBase()
94{
95     setType(12);
96     matrix_ = NULL;
97     startPositive_ = NULL;
98     startNegative_ = NULL;
99     lengths_ = NULL;
100     indices_ = NULL;
101     int iColumn;
102     assert (rhs.isColOrdered());
103     // get matrix data pointers
104     const int * row = rhs.getIndices();
105     const CoinBigIndex * columnStart = rhs.getVectorStarts();
106     const int * columnLength = rhs.getVectorLengths();
107     const double * elementByColumn = rhs.getElements();
108     numberColumns_ = rhs.getNumCols();
109     numberRows_ = -1;
110     indices_ = new int[rhs.getNumElements()];
111     startPositive_ = new CoinBigIndex [numberColumns_+1];
112     startNegative_ = new CoinBigIndex [numberColumns_];
113     int * temp = new int [rhs.getNumRows()];
114     CoinBigIndex j = 0;
115     CoinBigIndex numberGoodP = 0;
116     CoinBigIndex numberGoodM = 0;
117     CoinBigIndex numberBad = 0;
118     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
119          CoinBigIndex k;
120          int iNeg = 0;
121          startPositive_[iColumn] = j;
122          for (k = columnStart[iColumn]; k < columnStart[iColumn] + columnLength[iColumn];
123                    k++) {
124               int iRow;
125               if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
126                    iRow = row[k];
127                    numberRows_ = CoinMax(numberRows_, iRow);
128                    indices_[j++] = iRow;
129                    numberGoodP++;
130               } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
131                    iRow = row[k];
132                    numberRows_ = CoinMax(numberRows_, iRow);
133                    temp[iNeg++] = iRow;
134                    numberGoodM++;
135               } else {
136                    numberBad++;
137               }
138          }
139          // move negative
140          startNegative_[iColumn] = j;
141          for (k = 0; k < iNeg; k++) {
142               indices_[j++] = temp[k];
143          }
144     }
145     startPositive_[numberColumns_] = j;
146     delete [] temp;
147     if (numberBad) {
148          delete [] indices_;
149          indices_ = NULL;
150          numberRows_ = 0;
151          numberColumns_ = 0;
152          delete [] startPositive_;
153          delete [] startNegative_;
154          // Put in statistics
155          startPositive_ = new CoinBigIndex [3];
156          startPositive_[0] = numberGoodP;
157          startPositive_[1] = numberGoodM;
158          startPositive_[2] = numberBad;
159          startNegative_ = NULL;
160     } else {
161          numberRows_ ++; //  correct
162          // but number should be same as rhs
163          assert (numberRows_ <= rhs.getNumRows());
164          numberRows_ = rhs.getNumRows();
165          columnOrdered_ = true;
166     }
167     // Check valid
168     if (!numberBad)
169          checkValid(false);
170}
171
172//-------------------------------------------------------------------
173// Destructor
174//-------------------------------------------------------------------
175ClpPlusMinusOneMatrix::~ClpPlusMinusOneMatrix ()
176{
177     delete matrix_;
178     delete [] startPositive_;
179     delete [] startNegative_;
180     delete [] lengths_;
181     delete [] indices_;
182}
183
184//----------------------------------------------------------------
185// Assignment operator
186//-------------------------------------------------------------------
187ClpPlusMinusOneMatrix &
188ClpPlusMinusOneMatrix::operator=(const ClpPlusMinusOneMatrix& rhs)
189{
190     if (this != &rhs) {
191          ClpMatrixBase::operator=(rhs);
192          delete matrix_;
193          delete [] startPositive_;
194          delete [] startNegative_;
195          delete [] lengths_;
196          delete [] indices_;
197          matrix_ = NULL;
198          startPositive_ = NULL;
199          lengths_ = NULL;
200          indices_ = NULL;
201          numberRows_ = rhs.numberRows_;
202          numberColumns_ = rhs.numberColumns_;
203          columnOrdered_ = rhs.columnOrdered_;
204          if (numberColumns_) {
205               CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
206               indices_ = new int [ numberElements];
207               CoinMemcpyN(rhs.indices_, numberElements, indices_);
208               startPositive_ = new CoinBigIndex [ numberColumns_+1];
209               CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_);
210               startNegative_ = new CoinBigIndex [ numberColumns_];
211               CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_);
212          }
213     }
214     return *this;
215}
216//-------------------------------------------------------------------
217// Clone
218//-------------------------------------------------------------------
219ClpMatrixBase * ClpPlusMinusOneMatrix::clone() const
220{
221     return new ClpPlusMinusOneMatrix(*this);
222}
223/* Subset clone (without gaps).  Duplicates are allowed
224   and order is as given */
225ClpMatrixBase *
226ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
227                                    int numberColumns,
228                                    const int * whichColumns) const
229{
230     return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
231                                      numberColumns, whichColumns);
232}
233/* Subset constructor (without gaps).  Duplicates are allowed
234   and order is as given */
235ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
236     const ClpPlusMinusOneMatrix & rhs,
237     int numberRows, const int * whichRow,
238     int numberColumns, const int * whichColumn)
239     : ClpMatrixBase(rhs)
240{
241     matrix_ = NULL;
242     startPositive_ = NULL;
243     startNegative_ = NULL;
244     lengths_ = NULL;
245     indices_ = NULL;
246     numberRows_ = 0;
247     numberColumns_ = 0;
248     columnOrdered_ = rhs.columnOrdered_;
249     if (numberRows <= 0 || numberColumns <= 0) {
250          startPositive_ = new CoinBigIndex [1];
251          startPositive_[0] = 0;
252     } else {
253          numberColumns_ = numberColumns;
254          numberRows_ = numberRows;
255          const int * index1 = rhs.indices_;
256          CoinBigIndex * startPositive1 = rhs.startPositive_;
257
258          int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
259          int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
260          int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
261          int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
262          // Also swap incoming if not column ordered
263          if (!columnOrdered_) {
264               int temp1 = numberRows;
265               numberRows = numberColumns;
266               numberColumns = temp1;
267               const int * temp2;
268               temp2 = whichRow;
269               whichRow = whichColumn;
270               whichColumn = temp2;
271          }
272          // Throw exception if rhs empty
273          if (numberMajor1 <= 0 || numberMinor1 <= 0)
274               throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
275          // Array to say if an old row is in new copy
276          int * newRow = new int [numberMinor1];
277          int iRow;
278          for (iRow = 0; iRow < numberMinor1; iRow++)
279               newRow[iRow] = -1;
280          // and array for duplicating rows
281          int * duplicateRow = new int [numberMinor];
282          int numberBad = 0;
283          for (iRow = 0; iRow < numberMinor; iRow++) {
284               duplicateRow[iRow] = -1;
285               int kRow = whichRow[iRow];
286               if (kRow >= 0  && kRow < numberMinor1) {
287                    if (newRow[kRow] < 0) {
288                         // first time
289                         newRow[kRow] = iRow;
290                    } else {
291                         // duplicate
292                         int lastRow = newRow[kRow];
293                         newRow[kRow] = iRow;
294                         duplicateRow[iRow] = lastRow;
295                    }
296               } else {
297                    // bad row
298                    numberBad++;
299               }
300          }
301
302          if (numberBad)
303               throw CoinError("bad minor entries",
304                               "subset constructor", "ClpPlusMinusOneMatrix");
305          // now get size and check columns
306          CoinBigIndex size = 0;
307          int iColumn;
308          numberBad = 0;
309          for (iColumn = 0; iColumn < numberMajor; iColumn++) {
310               int kColumn = whichColumn[iColumn];
311               if (kColumn >= 0  && kColumn < numberMajor1) {
312                    CoinBigIndex i;
313                    for (i = startPositive1[kColumn]; i < startPositive1[kColumn+1]; i++) {
314                         int kRow = index1[i];
315                         kRow = newRow[kRow];
316                         while (kRow >= 0) {
317                              size++;
318                              kRow = duplicateRow[kRow];
319                         }
320                    }
321               } else {
322                    // bad column
323                    numberBad++;
324                    printf("%d %d %d %d\n", iColumn, numberMajor, numberMajor1, kColumn);
325               }
326          }
327          if (numberBad)
328               throw CoinError("bad major entries",
329                               "subset constructor", "ClpPlusMinusOneMatrix");
330          // now create arrays
331          startPositive_ = new CoinBigIndex [numberMajor+1];
332          startNegative_ = new CoinBigIndex [numberMajor];
333          indices_ = new int[size];
334          // and fill them
335          size = 0;
336          startPositive_[0] = 0;
337          CoinBigIndex * startNegative1 = rhs.startNegative_;
338          for (iColumn = 0; iColumn < numberMajor; iColumn++) {
339               int kColumn = whichColumn[iColumn];
340               CoinBigIndex i;
341               for (i = startPositive1[kColumn]; i < startNegative1[kColumn]; i++) {
342                    int kRow = index1[i];
343                    kRow = newRow[kRow];
344                    while (kRow >= 0) {
345                         indices_[size++] = kRow;
346                         kRow = duplicateRow[kRow];
347                    }
348               }
349               startNegative_[iColumn] = size;
350               for (; i < startPositive1[kColumn+1]; i++) {
351                    int kRow = index1[i];
352                    kRow = newRow[kRow];
353                    while (kRow >= 0) {
354                         indices_[size++] = kRow;
355                         kRow = duplicateRow[kRow];
356                    }
357               }
358               startPositive_[iColumn+1] = size;
359          }
360          delete [] newRow;
361          delete [] duplicateRow;
362     }
363     // Check valid
364     checkValid(false);
365}
366
367
368/* Returns a new matrix in reverse order without gaps */
369ClpMatrixBase *
370ClpPlusMinusOneMatrix::reverseOrderedCopy() const
371{
372     int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
373     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
374     // count number in each row/column
375     CoinBigIndex * tempP = new CoinBigIndex [numberMinor];
376     CoinBigIndex * tempN = new CoinBigIndex [numberMinor];
377     memset(tempP, 0, numberMinor * sizeof(CoinBigIndex));
378     memset(tempN, 0, numberMinor * sizeof(CoinBigIndex));
379     CoinBigIndex j = 0;
380     int i;
381     for (i = 0; i < numberMajor; i++) {
382          for (; j < startNegative_[i]; j++) {
383               int iRow = indices_[j];
384               tempP[iRow]++;
385          }
386          for (; j < startPositive_[i+1]; j++) {
387               int iRow = indices_[j];
388               tempN[iRow]++;
389          }
390     }
391     int * newIndices = new int [startPositive_[numberMajor]];
392     CoinBigIndex * newP = new CoinBigIndex [numberMinor+1];
393     CoinBigIndex * newN = new CoinBigIndex[numberMinor];
394     int iRow;
395     j = 0;
396     // do starts
397     for (iRow = 0; iRow < numberMinor; iRow++) {
398          newP[iRow] = j;
399          j += tempP[iRow];
400          tempP[iRow] = newP[iRow];
401          newN[iRow] = j;
402          j += tempN[iRow];
403          tempN[iRow] = newN[iRow];
404     }
405     newP[numberMinor] = j;
406     j = 0;
407     for (i = 0; i < numberMajor; i++) {
408          for (; j < startNegative_[i]; j++) {
409               int iRow = indices_[j];
410               CoinBigIndex put = tempP[iRow];
411               newIndices[put++] = i;
412               tempP[iRow] = put;
413          }
414          for (; j < startPositive_[i+1]; j++) {
415               int iRow = indices_[j];
416               CoinBigIndex put = tempN[iRow];
417               newIndices[put++] = i;
418               tempN[iRow] = put;
419          }
420     }
421     delete [] tempP;
422     delete [] tempN;
423     ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
424     newCopy->passInCopy(numberMinor, numberMajor,
425                         !columnOrdered_,  newIndices, newP, newN);
426     return newCopy;
427}
428//unscaled versions
429void
430ClpPlusMinusOneMatrix::times(double scalar,
431                             const double * x, double * y) const
432{
433     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
434     int i;
435     CoinBigIndex j;
436     assert (columnOrdered_);
437     for (i = 0; i < numberMajor; i++) {
438          double value = scalar * x[i];
439          if (value) {
440               for (j = startPositive_[i]; j < startNegative_[i]; j++) {
441                    int iRow = indices_[j];
442                    y[iRow] += value;
443               }
444               for (; j < startPositive_[i+1]; j++) {
445                    int iRow = indices_[j];
446                    y[iRow] -= value;
447               }
448          }
449     }
450}
451void
452ClpPlusMinusOneMatrix::transposeTimes(double scalar,
453                                      const double * x, double * y) const
454{
455     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
456     int i;
457     CoinBigIndex j = 0;
458     assert (columnOrdered_);
459     for (i = 0; i < numberMajor; i++) {
460          double value = 0.0;
461          for (; j < startNegative_[i]; j++) {
462               int iRow = indices_[j];
463               value += x[iRow];
464          }
465          for (; j < startPositive_[i+1]; j++) {
466               int iRow = indices_[j];
467               value -= x[iRow];
468          }
469          y[i] += scalar * value;
470     }
471}
472void
473ClpPlusMinusOneMatrix::times(double scalar,
474                             const double * x, double * y,
475                             const double * /*rowScale*/,
476                             const double * /*columnScale*/) const
477{
478     // we know it is not scaled
479     times(scalar, x, y);
480}
481void
482ClpPlusMinusOneMatrix::transposeTimes( double scalar,
483                                       const double * x, double * y,
484                                       const double * /*rowScale*/,
485                                       const double * /*columnScale*/,
486                                       double * /*spare*/) const
487{
488     // we know it is not scaled
489     transposeTimes(scalar, x, y);
490}
491/* Return <code>x * A + y</code> in <code>z</code>.
492        Squashes small elements and knows about ClpSimplex */
493void
494ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex * model, double scalar,
495                                      const CoinIndexedVector * rowArray,
496                                      CoinIndexedVector * y,
497                                      CoinIndexedVector * columnArray) const
498{
499     // we know it is not scaled
500     columnArray->clear();
501     double * pi = rowArray->denseVector();
502     int numberNonZero = 0;
503     int * index = columnArray->getIndices();
504     double * array = columnArray->denseVector();
505     int numberInRowArray = rowArray->getNumElements();
506     // maybe I need one in OsiSimplex
507     double zeroTolerance = model->zeroTolerance();
508     int numberRows = model->numberRows();
509     bool packed = rowArray->packedMode();
510#ifndef NO_RTTI
511     ClpPlusMinusOneMatrix* rowCopy =
512          dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
513#else
514     ClpPlusMinusOneMatrix* rowCopy =
515          static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
516#endif
517     double factor = 0.3;
518     // We may not want to do by row if there may be cache problems
519     int numberColumns = model->numberColumns();
520     // It would be nice to find L2 cache size - for moment 512K
521     // Be slightly optimistic
522     if (numberColumns * sizeof(double) > 1000000) {
523          if (numberRows * 10 < numberColumns)
524               factor = 0.1;
525          else if (numberRows * 4 < numberColumns)
526               factor = 0.15;
527          else if (numberRows * 2 < numberColumns)
528               factor = 0.2;
529     }
530     if (numberInRowArray > factor * numberRows || !rowCopy) {
531          assert (!y->getNumElements());
532          // do by column
533          // Need to expand if packed mode
534          int iColumn;
535          CoinBigIndex j = 0;
536          assert (columnOrdered_);
537          if (packed) {
538               // need to expand pi into y
539               assert(y->capacity() >= numberRows);
540               double * piOld = pi;
541               pi = y->denseVector();
542               const int * whichRow = rowArray->getIndices();
543               int i;
544               // modify pi so can collapse to one loop
545               for (i = 0; i < numberInRowArray; i++) {
546                    int iRow = whichRow[i];
547                    pi[iRow] = scalar * piOld[i];
548               }
549               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
550                    double value = 0.0;
551                    for (; j < startNegative_[iColumn]; j++) {
552                         int iRow = indices_[j];
553                         value += pi[iRow];
554                    }
555                    for (; j < startPositive_[iColumn+1]; j++) {
556                         int iRow = indices_[j];
557                         value -= pi[iRow];
558                    }
559                    if (fabs(value) > zeroTolerance) {
560                         array[numberNonZero] = value;
561                         index[numberNonZero++] = iColumn;
562                    }
563               }
564               for (i = 0; i < numberInRowArray; i++) {
565                    int iRow = whichRow[i];
566                    pi[iRow] = 0.0;
567               }
568          } else {
569               for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
570                    double value = 0.0;
571                    for (; j < startNegative_[iColumn]; j++) {
572                         int iRow = indices_[j];
573                         value += pi[iRow];
574                    }
575                    for (; j < startPositive_[iColumn+1]; j++) {
576                         int iRow = indices_[j];
577                         value -= pi[iRow];
578                    }
579                    value *= scalar;
580                    if (fabs(value) > zeroTolerance) {
581                         index[numberNonZero++] = iColumn;
582                         array[iColumn] = value;
583                    }
584               }
585          }
586          columnArray->setNumElements(numberNonZero);
587     } else {
588          // do by row
589          rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
590     }
591}
592/* Return <code>x * A + y</code> in <code>z</code>.
593        Squashes small elements and knows about ClpSimplex */
594void
595ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
596          const CoinIndexedVector * rowArray,
597          CoinIndexedVector * y,
598          CoinIndexedVector * columnArray) const
599{
600     columnArray->clear();
601     double * pi = rowArray->denseVector();
602     int numberNonZero = 0;
603     int * index = columnArray->getIndices();
604     double * array = columnArray->denseVector();
605     int numberInRowArray = rowArray->getNumElements();
606     // maybe I need one in OsiSimplex
607     double zeroTolerance = model->zeroTolerance();
608     const int * column = indices_;
609     const CoinBigIndex * startPositive = startPositive_;
610     const CoinBigIndex * startNegative = startNegative_;
611     const int * whichRow = rowArray->getIndices();
612     bool packed = rowArray->packedMode();
613     if (numberInRowArray > 2) {
614          // do by rows
615          int iRow;
616          double * markVector = y->denseVector(); // probably empty .. but
617          int numberOriginal = 0;
618          int i;
619          if (packed) {
620               numberNonZero = 0;
621               // and set up mark as char array
622               char * marked = reinterpret_cast<char *> (index + columnArray->capacity());
623               double * array2 = y->denseVector();
624#ifdef CLP_DEBUG
625               int numberColumns = model->numberColumns();
626               for (i = 0; i < numberColumns; i++) {
627                    assert(!marked[i]);
628                    assert(!array2[i]);
629               }
630#endif
631               for (i = 0; i < numberInRowArray; i++) {
632                    iRow = whichRow[i];
633                    double value = pi[i] * scalar;
634                    CoinBigIndex j;
635                    for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
636                         int iColumn = column[j];
637                         if (!marked[iColumn]) {
638                              marked[iColumn] = 1;
639                              index[numberNonZero++] = iColumn;
640                         }
641                         array2[iColumn] += value;
642                    }
643                    for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
644                         int iColumn = column[j];
645                         if (!marked[iColumn]) {
646                              marked[iColumn] = 1;
647                              index[numberNonZero++] = iColumn;
648                         }
649                         array2[iColumn] -= value;
650                    }
651               }
652               // get rid of tiny values and zero out marked
653               numberOriginal = numberNonZero;
654               numberNonZero = 0;
655               for (i = 0; i < numberOriginal; i++) {
656                    int iColumn = index[i];
657                    if (marked[iColumn]) {
658                         double value = array2[iColumn];
659                         array2[iColumn] = 0.0;
660                         marked[iColumn] = 0;
661                         if (fabs(value) > zeroTolerance) {
662                              array[numberNonZero] = value;
663                              index[numberNonZero++] = iColumn;
664                         }
665                    }
666               }
667          } else {
668               numberNonZero = 0;
669               // and set up mark as char array
670               char * marked = reinterpret_cast<char *> (markVector);
671               for (i = 0; i < numberOriginal; i++) {
672                    int iColumn = index[i];
673                    marked[iColumn] = 0;
674               }
675               for (i = 0; i < numberInRowArray; i++) {
676                    iRow = whichRow[i];
677                    double value = pi[iRow] * scalar;
678                    CoinBigIndex j;
679                    for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
680                         int iColumn = column[j];
681                         if (!marked[iColumn]) {
682                              marked[iColumn] = 1;
683                              index[numberNonZero++] = iColumn;
684                         }
685                         array[iColumn] += value;
686                    }
687                    for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
688                         int iColumn = column[j];
689                         if (!marked[iColumn]) {
690                              marked[iColumn] = 1;
691                              index[numberNonZero++] = iColumn;
692                         }
693                         array[iColumn] -= value;
694                    }
695               }
696               // get rid of tiny values and zero out marked
697               numberOriginal = numberNonZero;
698               numberNonZero = 0;
699               for (i = 0; i < numberOriginal; i++) {
700                    int iColumn = index[i];
701                    marked[iColumn] = 0;
702                    if (fabs(array[iColumn]) > zeroTolerance) {
703                         index[numberNonZero++] = iColumn;
704                    } else {
705                         array[iColumn] = 0.0;
706                    }
707               }
708          }
709     } else if (numberInRowArray == 2) {
710          /* do by rows when two rows (do longer first when not packed
711             and shorter first if packed */
712          int iRow0 = whichRow[0];
713          int iRow1 = whichRow[1];
714          CoinBigIndex j;
715          if (packed) {
716               double pi0 = pi[0];
717               double pi1 = pi[1];
718               if (startPositive[iRow0+1] - startPositive[iRow0] >
719                         startPositive[iRow1+1] - startPositive[iRow1]) {
720                    int temp = iRow0;
721                    iRow0 = iRow1;
722                    iRow1 = temp;
723                    pi0 = pi1;
724                    pi1 = pi[0];
725               }
726               // and set up mark as char array
727               char * marked = reinterpret_cast<char *> (index + columnArray->capacity());
728               int * lookup = y->getIndices();
729               double value = pi0 * scalar;
730               for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
731                    int iColumn = column[j];
732                    array[numberNonZero] = value;
733                    marked[iColumn] = 1;
734                    lookup[iColumn] = numberNonZero;
735                    index[numberNonZero++] = iColumn;
736               }
737               for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
738                    int iColumn = column[j];
739                    array[numberNonZero] = -value;
740                    marked[iColumn] = 1;
741                    lookup[iColumn] = numberNonZero;
742                    index[numberNonZero++] = iColumn;
743               }
744               int numberOriginal = numberNonZero;
745               value = pi1 * scalar;
746               for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
747                    int iColumn = column[j];
748                    if (marked[iColumn]) {
749                         int iLookup = lookup[iColumn];
750                         array[iLookup] += value;
751                    } else {
752                         if (fabs(value) > zeroTolerance) {
753                              array[numberNonZero] = value;
754                              index[numberNonZero++] = iColumn;
755                         }
756                    }
757               }
758               for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
759                    int iColumn = column[j];
760                    if (marked[iColumn]) {
761                         int iLookup = lookup[iColumn];
762                         array[iLookup] -= value;
763                    } else {
764                         if (fabs(value) > zeroTolerance) {
765                              array[numberNonZero] = -value;
766                              index[numberNonZero++] = iColumn;
767                         }
768                    }
769               }
770               // get rid of tiny values and zero out marked
771               int nDelete = 0;
772               for (j = 0; j < numberOriginal; j++) {
773                    int iColumn = index[j];
774                    marked[iColumn] = 0;
775                    if (fabs(array[j]) <= zeroTolerance)
776                         nDelete++;
777               }
778               if (nDelete) {
779                    numberOriginal = numberNonZero;
780                    numberNonZero = 0;
781                    for (j = 0; j < numberOriginal; j++) {
782                         int iColumn = index[j];
783                         double value = array[j];
784                         array[j] = 0.0;
785                         if (fabs(value) > zeroTolerance) {
786                              array[numberNonZero] = value;
787                              index[numberNonZero++] = iColumn;
788                         }
789                    }
790               }
791          } else {
792               if (startPositive[iRow0+1] - startPositive[iRow0] <
793                         startPositive[iRow1+1] - startPositive[iRow1]) {
794                    int temp = iRow0;
795                    iRow0 = iRow1;
796                    iRow1 = temp;
797               }
798               int numberOriginal;
799               int i;
800               numberNonZero = 0;
801               double value;
802               value = pi[iRow0] * scalar;
803               CoinBigIndex j;
804               for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
805                    int iColumn = column[j];
806                    index[numberNonZero++] = iColumn;
807                    array[iColumn] = value;
808               }
809               for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
810                    int iColumn = column[j];
811                    index[numberNonZero++] = iColumn;
812                    array[iColumn] = -value;
813               }
814               value = pi[iRow1] * scalar;
815               for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
816                    int iColumn = column[j];
817                    double value2 = array[iColumn];
818                    if (value2) {
819                         value2 += value;
820                    } else {
821                         value2 = value;
822                         index[numberNonZero++] = iColumn;
823                    }
824                    array[iColumn] = value2;
825               }
826               for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
827                    int iColumn = column[j];
828                    double value2 = array[iColumn];
829                    if (value2) {
830                         value2 -= value;
831                    } else {
832                         value2 = -value;
833                         index[numberNonZero++] = iColumn;
834                    }
835                    array[iColumn] = value2;
836               }
837               // get rid of tiny values and zero out marked
838               numberOriginal = numberNonZero;
839               numberNonZero = 0;
840               for (i = 0; i < numberOriginal; i++) {
841                    int iColumn = index[i];
842                    if (fabs(array[iColumn]) > zeroTolerance) {
843                         index[numberNonZero++] = iColumn;
844                    } else {
845                         array[iColumn] = 0.0;
846                    }
847               }
848          }
849     } else if (numberInRowArray == 1) {
850          // Just one row
851          int iRow = rowArray->getIndices()[0];
852          numberNonZero = 0;
853          double value;
854          iRow = whichRow[0];
855          CoinBigIndex j;
856          if (packed) {
857               value = pi[0] * scalar;
858               if (fabs(value) > zeroTolerance) {
859                    for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
860                         int iColumn = column[j];
861                         array[numberNonZero] = value;
862                         index[numberNonZero++] = iColumn;
863                    }
864                    for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
865                         int iColumn = column[j];
866                         array[numberNonZero] = -value;
867                         index[numberNonZero++] = iColumn;
868                    }
869               }
870          } else {
871               value = pi[iRow] * scalar;
872               if (fabs(value) > zeroTolerance) {
873                    for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
874                         int iColumn = column[j];
875                         array[iColumn] = value;
876                         index[numberNonZero++] = iColumn;
877                    }
878                    for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
879                         int iColumn = column[j];
880                         array[iColumn] = -value;
881                         index[numberNonZero++] = iColumn;
882                    }
883               }
884          }
885     }
886     columnArray->setNumElements(numberNonZero);
887     if (packed)
888          columnArray->setPacked();
889     y->setNumElements(0);
890}
891/* Return <code>x *A in <code>z</code> but
892   just for indices in y. */
893void
894ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex * ,
895          const CoinIndexedVector * rowArray,
896          const CoinIndexedVector * y,
897          CoinIndexedVector * columnArray) const
898{
899     columnArray->clear();
900     double * pi = rowArray->denseVector();
901     double * array = columnArray->denseVector();
902     int jColumn;
903     int numberToDo = y->getNumElements();
904     const int * which = y->getIndices();
905     assert (!rowArray->packedMode());
906     columnArray->setPacked();
907     for (jColumn = 0; jColumn < numberToDo; jColumn++) {
908          int iColumn = which[jColumn];
909          double value = 0.0;
910          CoinBigIndex j = startPositive_[iColumn];
911          for (; j < startNegative_[iColumn]; j++) {
912               int iRow = indices_[j];
913               value += pi[iRow];
914          }
915          for (; j < startPositive_[iColumn+1]; j++) {
916               int iRow = indices_[j];
917               value -= pi[iRow];
918          }
919          array[jColumn] = value;
920     }
921}
922/// returns number of elements in column part of basis,
923CoinBigIndex
924ClpPlusMinusOneMatrix::countBasis(const int * whichColumn,
925                                  int & numberColumnBasic)
926{
927     int i;
928     CoinBigIndex numberElements = 0;
929     for (i = 0; i < numberColumnBasic; i++) {
930          int iColumn = whichColumn[i];
931          numberElements += startPositive_[iColumn+1] - startPositive_[iColumn];
932     }
933     return numberElements;
934}
935void
936ClpPlusMinusOneMatrix::fillBasis(ClpSimplex * ,
937                                 const int * whichColumn,
938                                 int & numberColumnBasic,
939                                 int * indexRowU, int * start,
940                                 int * rowCount, int * columnCount,
941                                 CoinFactorizationDouble * elementU)
942{
943     int i;
944     CoinBigIndex numberElements = start[0];
945     assert (columnOrdered_);
946     for (i = 0; i < numberColumnBasic; i++) {
947          int iColumn = whichColumn[i];
948          CoinBigIndex j = startPositive_[iColumn];
949          for (; j < startNegative_[iColumn]; j++) {
950               int iRow = indices_[j];
951               indexRowU[numberElements] = iRow;
952               rowCount[iRow]++;
953               elementU[numberElements++] = 1.0;
954          }
955          for (; j < startPositive_[iColumn+1]; j++) {
956               int iRow = indices_[j];
957               indexRowU[numberElements] = iRow;
958               rowCount[iRow]++;
959               elementU[numberElements++] = -1.0;
960          }
961          start[i+1] = numberElements;
962          columnCount[i] = numberElements - start[i];
963     }
964}
965/* Unpacks a column into an CoinIndexedvector
966 */
967void
968ClpPlusMinusOneMatrix::unpack(const ClpSimplex * ,
969                              CoinIndexedVector * rowArray,
970                              int iColumn) const
971{
972     CoinBigIndex j = startPositive_[iColumn];
973     for (; j < startNegative_[iColumn]; j++) {
974          int iRow = indices_[j];
975          rowArray->add(iRow, 1.0);
976     }
977     for (; j < startPositive_[iColumn+1]; j++) {
978          int iRow = indices_[j];
979          rowArray->add(iRow, -1.0);
980     }
981}
982/* Unpacks a column into an CoinIndexedvector
983** in packed foramt
984Note that model is NOT const.  Bounds and objective could
985be modified if doing column generation (just for this variable) */
986void
987ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * ,
988                                    CoinIndexedVector * rowArray,
989                                    int iColumn) const
990{
991     int * index = rowArray->getIndices();
992     double * array = rowArray->denseVector();
993     int number = 0;
994     CoinBigIndex j = startPositive_[iColumn];
995     for (; j < startNegative_[iColumn]; j++) {
996          int iRow = indices_[j];
997          array[number] = 1.0;
998          index[number++] = iRow;
999     }
1000     for (; j < startPositive_[iColumn+1]; j++) {
1001          int iRow = indices_[j];
1002          array[number] = -1.0;
1003          index[number++] = iRow;
1004     }
1005     rowArray->setNumElements(number);
1006     rowArray->setPackedMode(true);
1007}
1008/* Adds multiple of a column into an CoinIndexedvector
1009      You can use quickAdd to add to vector */
1010void
1011ClpPlusMinusOneMatrix::add(const ClpSimplex * , CoinIndexedVector * rowArray,
1012                           int iColumn, double multiplier) const
1013{
1014     CoinBigIndex j = startPositive_[iColumn];
1015     for (; j < startNegative_[iColumn]; j++) {
1016          int iRow = indices_[j];
1017          rowArray->quickAdd(iRow, multiplier);
1018     }
1019     for (; j < startPositive_[iColumn+1]; j++) {
1020          int iRow = indices_[j];
1021          rowArray->quickAdd(iRow, -multiplier);
1022     }
1023}
1024/* Adds multiple of a column into an array */
1025void
1026ClpPlusMinusOneMatrix::add(const ClpSimplex * , double * array,
1027                           int iColumn, double multiplier) const
1028{
1029     CoinBigIndex j = startPositive_[iColumn];
1030     for (; j < startNegative_[iColumn]; j++) {
1031          int iRow = indices_[j];
1032          array[iRow] += multiplier;
1033     }
1034     for (; j < startPositive_[iColumn+1]; j++) {
1035          int iRow = indices_[j];
1036          array[iRow] -= multiplier;
1037     }
1038}
1039
1040// Return a complete CoinPackedMatrix
1041CoinPackedMatrix *
1042ClpPlusMinusOneMatrix::getPackedMatrix() const
1043{
1044     if (!matrix_) {
1045          int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
1046          int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1047          int numberElements = startPositive_[numberMajor];
1048          double * elements = new double [numberElements];
1049          CoinBigIndex j = 0;
1050          int i;
1051          for (i = 0; i < numberMajor; i++) {
1052               for (; j < startNegative_[i]; j++) {
1053                    elements[j] = 1.0;
1054               }
1055               for (; j < startPositive_[i+1]; j++) {
1056                    elements[j] = -1.0;
1057               }
1058          }
1059          matrix_ =  new CoinPackedMatrix(columnOrdered_, numberMinor, numberMajor,
1060                                          getNumElements(),
1061                                          elements, indices_,
1062                                          startPositive_, getVectorLengths());
1063          delete [] elements;
1064          delete [] lengths_;
1065          lengths_ = NULL;
1066     }
1067     return matrix_;
1068}
1069/* A vector containing the elements in the packed matrix. Note that there
1070   might be gaps in this list, entries that do not belong to any
1071   major-dimension vector. To get the actual elements one should look at
1072   this vector together with vectorStarts and vectorLengths. */
1073const double *
1074ClpPlusMinusOneMatrix::getElements() const
1075{
1076     if (!matrix_)
1077          getPackedMatrix();
1078     return matrix_->getElements();
1079}
1080
1081const CoinBigIndex *
1082ClpPlusMinusOneMatrix::getVectorStarts() const
1083{
1084     return startPositive_;
1085}
1086/* The lengths of the major-dimension vectors. */
1087const int *
1088ClpPlusMinusOneMatrix::getVectorLengths() const
1089{
1090     if (!lengths_) {
1091          int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1092          lengths_ = new int [numberMajor];
1093          int i;
1094          for (i = 0; i < numberMajor; i++) {
1095               lengths_[i] = startPositive_[i+1] - startPositive_[i];
1096          }
1097     }
1098     return lengths_;
1099}
1100/* Delete the columns whose indices are listed in <code>indDel</code>. */
1101void
1102ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel)
1103{
1104     int iColumn;
1105     CoinBigIndex newSize = startPositive_[numberColumns_];;
1106     int numberBad = 0;
1107     // Use array to make sure we can have duplicates
1108     int * which = new int[numberColumns_];
1109     memset(which, 0, numberColumns_ * sizeof(int));
1110     int nDuplicate = 0;
1111     for (iColumn = 0; iColumn < numDel; iColumn++) {
1112          int jColumn = indDel[iColumn];
1113          if (jColumn < 0 || jColumn >= numberColumns_) {
1114               numberBad++;
1115          } else {
1116               newSize -= startPositive_[jColumn+1] - startPositive_[jColumn];
1117               if (which[jColumn])
1118                    nDuplicate++;
1119               else
1120                    which[jColumn] = 1;
1121          }
1122     }
1123     if (numberBad)
1124          throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
1125     int newNumber = numberColumns_ - numDel + nDuplicate;
1126     // Get rid of temporary arrays
1127     delete [] lengths_;
1128     lengths_ = NULL;
1129     delete matrix_;
1130     matrix_ = NULL;
1131     CoinBigIndex * newPositive = new CoinBigIndex [newNumber+1];
1132     CoinBigIndex * newNegative = new CoinBigIndex [newNumber];
1133     int * newIndices = new int [newSize];
1134     newNumber = 0;
1135     newSize = 0;
1136     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1137          if (!which[iColumn]) {
1138               CoinBigIndex start, end;
1139               CoinBigIndex i;
1140               start = startPositive_[iColumn];
1141               end = startNegative_[iColumn];
1142               newPositive[newNumber] = newSize;
1143               for (i = start; i < end; i++)
1144                    newIndices[newSize++] = indices_[i];
1145               start = startNegative_[iColumn];
1146               end = startPositive_[iColumn+1];
1147               newNegative[newNumber++] = newSize;
1148               for (i = start; i < end; i++)
1149                    newIndices[newSize++] = indices_[i];
1150          }
1151     }
1152     newPositive[newNumber] = newSize;
1153     delete [] which;
1154     delete [] startPositive_;
1155     startPositive_ = newPositive;
1156     delete [] startNegative_;
1157     startNegative_ = newNegative;
1158     delete [] indices_;
1159     indices_ = newIndices;
1160     numberColumns_ = newNumber;
1161}
1162/* Delete the rows whose indices are listed in <code>indDel</code>. */
1163void
1164ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel)
1165{
1166     int iRow;
1167     int numberBad = 0;
1168     // Use array to make sure we can have duplicates
1169     int * which = new int[numberRows_];
1170     memset(which, 0, numberRows_ * sizeof(int));
1171     int nDuplicate = 0;
1172     for (iRow = 0; iRow < numDel; iRow++) {
1173          int jRow = indDel[iRow];
1174          if (jRow < 0 || jRow >= numberRows_) {
1175               numberBad++;
1176          } else {
1177               if (which[jRow])
1178                    nDuplicate++;
1179               else
1180                    which[jRow] = 1;
1181          }
1182     }
1183     if (numberBad)
1184          throw CoinError("Indices out of range", "deleteRows", "ClpPlusMinusOneMatrix");
1185     CoinBigIndex iElement;
1186     CoinBigIndex numberElements = startPositive_[numberColumns_];
1187     CoinBigIndex newSize = 0;
1188     for (iElement = 0; iElement < numberElements; iElement++) {
1189          iRow = indices_[iElement];
1190          if (!which[iRow])
1191               newSize++;
1192     }
1193     int newNumber = numberRows_ - numDel + nDuplicate;
1194     // Get rid of temporary arrays
1195     delete [] lengths_;
1196     lengths_ = NULL;
1197     delete matrix_;
1198     matrix_ = NULL;
1199     int * newIndices = new int [newSize];
1200     newSize = 0;
1201     int iColumn;
1202     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1203          CoinBigIndex start, end;
1204          CoinBigIndex i;
1205          start = startPositive_[iColumn];
1206          end = startNegative_[iColumn];
1207          startPositive_[newNumber] = newSize;
1208          for (i = start; i < end; i++) {
1209               iRow = indices_[i];
1210               if (!which[iRow])
1211                    newIndices[newSize++] = iRow;
1212          }
1213          start = startNegative_[iColumn];
1214          end = startPositive_[iColumn+1];
1215          startNegative_[newNumber] = newSize;
1216          for (i = start; i < end; i++) {
1217               iRow = indices_[i];
1218               if (!which[iRow])
1219                    newIndices[newSize++] = iRow;
1220          }
1221     }
1222     startPositive_[numberColumns_] = newSize;
1223     delete [] which;
1224     delete [] indices_;
1225     indices_ = newIndices;
1226     numberRows_ = newNumber;
1227}
1228bool
1229ClpPlusMinusOneMatrix::isColOrdered() const
1230{
1231     return columnOrdered_;
1232}
1233/* Number of entries in the packed matrix. */
1234CoinBigIndex
1235ClpPlusMinusOneMatrix::getNumElements() const
1236{
1237     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1238     if (startPositive_)
1239          return startPositive_[numberMajor];
1240     else
1241          return 0;
1242}
1243// pass in copy (object takes ownership)
1244void
1245ClpPlusMinusOneMatrix::passInCopy(int numberRows, int numberColumns,
1246                                  bool columnOrdered, int * indices,
1247                                  CoinBigIndex * startPositive, CoinBigIndex * startNegative)
1248{
1249     columnOrdered_ = columnOrdered;
1250     startPositive_ = startPositive;
1251     startNegative_ = startNegative;
1252     indices_ = indices;
1253     numberRows_ = numberRows;
1254     numberColumns_ = numberColumns;
1255     // Check valid
1256     checkValid(false);
1257}
1258// Just checks matrix valid - will say if dimensions not quite right if detail
1259void
1260ClpPlusMinusOneMatrix::checkValid(bool detail) const
1261{
1262     int maxIndex = -1;
1263     int minIndex = columnOrdered_ ? numberRows_ : numberColumns_;
1264     int number = !columnOrdered_ ? numberRows_ : numberColumns_;
1265     int numberElements = getNumElements();
1266     CoinBigIndex last = -1;
1267     int bad = 0;
1268     for (int i = 0; i < number; i++) {
1269          if(startPositive_[i] < last)
1270               bad++;
1271          else
1272               last = startPositive_[i];
1273          if(startNegative_[i] < last)
1274               bad++;
1275          else
1276               last = startNegative_[i];
1277     }
1278     if(startPositive_[number] < last)
1279          bad++;
1280     CoinAssertHint(!bad, "starts are not monotonic");
1281     for (CoinBigIndex cbi = 0; cbi < numberElements; cbi++) {
1282          maxIndex = CoinMax(indices_[cbi], maxIndex);
1283          minIndex = CoinMin(indices_[cbi], minIndex);
1284     }
1285     CoinAssert(maxIndex < (columnOrdered_ ? numberRows_ : numberColumns_));
1286     CoinAssert(minIndex >= 0);
1287     if (detail) {
1288          if (minIndex > 0 || maxIndex + 1 < (columnOrdered_ ? numberRows_ : numberColumns_))
1289               printf("Not full range of indices - %d to %d\n", minIndex, maxIndex);
1290     }
1291}
1292/* Given positive integer weights for each row fills in sum of weights
1293   for each column (and slack).
1294   Returns weights vector
1295*/
1296CoinBigIndex *
1297ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
1298{
1299     int numberRows = model->numberRows();
1300     int numberColumns = model->numberColumns();
1301     int number = numberRows + numberColumns;
1302     CoinBigIndex * weights = new CoinBigIndex[number];
1303     int i;
1304     for (i = 0; i < numberColumns; i++) {
1305          CoinBigIndex j;
1306          CoinBigIndex count = 0;
1307          for (j = startPositive_[i]; j < startPositive_[i+1]; j++) {
1308               int iRow = indices_[j];
1309               count += inputWeights[iRow];
1310          }
1311          weights[i] = count;
1312     }
1313     for (i = 0; i < numberRows; i++) {
1314          weights[i+numberColumns] = inputWeights[i];
1315     }
1316     return weights;
1317}
1318// Append Columns
1319void
1320ClpPlusMinusOneMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
1321{
1322     int iColumn;
1323     CoinBigIndex size = 0;
1324     int numberBad = 0;
1325     for (iColumn = 0; iColumn < number; iColumn++) {
1326          int n = columns[iColumn]->getNumElements();
1327          const double * element = columns[iColumn]->getElements();
1328          size += n;
1329          int i;
1330          for (i = 0; i < n; i++) {
1331               if (fabs(element[i]) != 1.0)
1332                    numberBad++;
1333          }
1334     }
1335     if (numberBad)
1336          throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix");
1337     // Get rid of temporary arrays
1338     delete [] lengths_;
1339     lengths_ = NULL;
1340     delete matrix_;
1341     matrix_ = NULL;
1342     int numberNow = startPositive_[numberColumns_];
1343     CoinBigIndex * temp;
1344     temp = new CoinBigIndex [numberColumns_+1+number];
1345     CoinMemcpyN(startPositive_, (numberColumns_ + 1), temp);
1346     delete [] startPositive_;
1347     startPositive_ = temp;
1348     temp = new CoinBigIndex [numberColumns_+number];
1349     CoinMemcpyN(startNegative_, numberColumns_, temp);
1350     delete [] startNegative_;
1351     startNegative_ = temp;
1352     int * temp2 = new int [numberNow+size];
1353     CoinMemcpyN(indices_, numberNow, temp2);
1354     delete [] indices_;
1355     indices_ = temp2;
1356     // now add
1357     size = numberNow;
1358     for (iColumn = 0; iColumn < number; iColumn++) {
1359          int n = columns[iColumn]->getNumElements();
1360          const int * row = columns[iColumn]->getIndices();
1361          const double * element = columns[iColumn]->getElements();
1362          int i;
1363          for (i = 0; i < n; i++) {
1364               if (element[i] == 1.0)
1365                    indices_[size++] = row[i];
1366          }
1367          startNegative_[iColumn+numberColumns_] = size;
1368          for (i = 0; i < n; i++) {
1369               if (element[i] == -1.0)
1370                    indices_[size++] = row[i];
1371          }
1372          startPositive_[iColumn+numberColumns_+1] = size;
1373     }
1374
1375     numberColumns_ += number;
1376}
1377// Append Rows
1378void
1379ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
1380{
1381     // Allocate arrays to use for counting
1382     int * countPositive = new int [numberColumns_+1];
1383     memset(countPositive, 0, numberColumns_ * sizeof(int));
1384     int * countNegative = new int [numberColumns_];
1385     memset(countNegative, 0, numberColumns_ * sizeof(int));
1386     int iRow;
1387     CoinBigIndex size = 0;
1388     int numberBad = 0;
1389     for (iRow = 0; iRow < number; iRow++) {
1390          int n = rows[iRow]->getNumElements();
1391          const int * column = rows[iRow]->getIndices();
1392          const double * element = rows[iRow]->getElements();
1393          size += n;
1394          int i;
1395          for (i = 0; i < n; i++) {
1396               int iColumn = column[i];
1397               if (element[i] == 1.0)
1398                    countPositive[iColumn]++;
1399               else if (element[i] == -1.0)
1400                    countNegative[iColumn]++;
1401               else
1402                    numberBad++;
1403          }
1404     }
1405     if (numberBad)
1406          throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix");
1407     // Get rid of temporary arrays
1408     delete [] lengths_;
1409     lengths_ = NULL;
1410     delete matrix_;
1411     matrix_ = NULL;
1412     int numberNow = startPositive_[numberColumns_];
1413     int * newIndices = new int [numberNow+size];
1414     // Update starts and turn counts into positions
1415     // also move current indices
1416     int iColumn;
1417     CoinBigIndex numberAdded = 0;
1418     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1419          int n, move;
1420          CoinBigIndex now;
1421          now = startPositive_[iColumn];
1422          move = startNegative_[iColumn] - now;
1423          n = countPositive[iColumn];
1424          startPositive_[iColumn] += numberAdded;
1425          CoinMemcpyN(newIndices + startPositive_[iColumn], move, indices_ + now);
1426          countPositive[iColumn] = startNegative_[iColumn] + numberAdded;
1427          numberAdded += n;
1428          now = startNegative_[iColumn];
1429          move = startPositive_[iColumn+1] - now;
1430          n = countNegative[iColumn];
1431          startNegative_[iColumn] += numberAdded;
1432          CoinMemcpyN(newIndices + startNegative_[iColumn], move, indices_ + now);
1433          countNegative[iColumn] = startPositive_[iColumn+1] + numberAdded;
1434          numberAdded += n;
1435     }
1436     delete [] indices_;
1437     indices_ = newIndices;
1438     startPositive_[numberColumns_] += numberAdded;
1439     // Now put in
1440     for (iRow = 0; iRow < number; iRow++) {
1441          int newRow = numberRows_ + iRow;
1442          int n = rows[iRow]->getNumElements();
1443          const int * column = rows[iRow]->getIndices();
1444          const double * element = rows[iRow]->getElements();
1445          int i;
1446          for (i = 0; i < n; i++) {
1447               int iColumn = column[i];
1448               int put;
1449               if (element[i] == 1.0) {
1450                    put = countPositive[iColumn];
1451                    countPositive[iColumn] = put + 1;
1452               } else {
1453                    put = countNegative[iColumn];
1454                    countNegative[iColumn] = put + 1;
1455               }
1456               indices_[put] = newRow;
1457          }
1458     }
1459     delete [] countPositive;
1460     delete [] countNegative;
1461     numberRows_ += number;
1462}
1463/* Returns largest and smallest elements of both signs.
1464   Largest refers to largest absolute value.
1465*/
1466void
1467ClpPlusMinusOneMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
1468                                       double & smallestPositive, double & largestPositive)
1469{
1470     int iColumn;
1471     bool plusOne = false;
1472     bool minusOne = false;
1473     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1474          if (startNegative_[iColumn] > startPositive_[iColumn])
1475               plusOne = true;
1476          if (startPositive_[iColumn+1] > startNegative_[iColumn])
1477               minusOne = true;
1478     }
1479     if (minusOne) {
1480          smallestNegative = -1.0;
1481          largestNegative = -1.0;
1482     } else {
1483          smallestNegative = 0.0;
1484          largestNegative = 0.0;
1485     }
1486     if (plusOne) {
1487          smallestPositive = 1.0;
1488          largestPositive = 1.0;
1489     } else {
1490          smallestPositive = 0.0;
1491          largestPositive = 0.0;
1492     }
1493}
1494// Says whether it can do partial pricing
1495bool
1496ClpPlusMinusOneMatrix::canDoPartialPricing() const
1497{
1498     return true;
1499}
1500// Partial pricing
1501void
1502ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
1503                                      int & bestSequence, int & numberWanted)
1504{
1505     numberWanted = currentWanted_;
1506     int start = static_cast<int> (startFraction * numberColumns_);
1507     int end = CoinMin(static_cast<int> (endFraction * numberColumns_ + 1), numberColumns_);
1508     CoinBigIndex j;
1509     double tolerance = model->currentDualTolerance();
1510     double * reducedCost = model->djRegion();
1511     const double * duals = model->dualRowSolution();
1512     const double * cost = model->costRegion();
1513     double bestDj;
1514     if (bestSequence >= 0)
1515          bestDj = fabs(reducedCost[bestSequence]);
1516     else
1517          bestDj = tolerance;
1518     int sequenceOut = model->sequenceOut();
1519     int saveSequence = bestSequence;
1520     int iSequence;
1521     for (iSequence = start; iSequence < end; iSequence++) {
1522          if (iSequence != sequenceOut) {
1523               double value;
1524               ClpSimplex::Status status = model->getStatus(iSequence);
1525
1526               switch(status) {
1527
1528               case ClpSimplex::basic:
1529               case ClpSimplex::isFixed:
1530                    break;
1531               case ClpSimplex::isFree:
1532               case ClpSimplex::superBasic:
1533                    value = cost[iSequence];
1534                    j = startPositive_[iSequence];
1535                    for (; j < startNegative_[iSequence]; j++) {
1536                         int iRow = indices_[j];
1537                         value -= duals[iRow];
1538                    }
1539                    for (; j < startPositive_[iSequence+1]; j++) {
1540                         int iRow = indices_[j];
1541                         value += duals[iRow];
1542                    }
1543                    value = fabs(value);
1544                    if (value > FREE_ACCEPT * tolerance) {
1545                         numberWanted--;
1546                         // we are going to bias towards free (but only if reasonable)
1547                         value *= FREE_BIAS;
1548                         if (value > bestDj) {
1549                              // check flagged variable and correct dj
1550                              if (!model->flagged(iSequence)) {
1551                                   bestDj = value;
1552                                   bestSequence = iSequence;
1553                              } else {
1554                                   // just to make sure we don't exit before got something
1555                                   numberWanted++;
1556                              }
1557                         }
1558                    }
1559                    break;
1560               case ClpSimplex::atUpperBound:
1561                    value = cost[iSequence];
1562                    j = startPositive_[iSequence];
1563                    for (; j < startNegative_[iSequence]; j++) {
1564                         int iRow = indices_[j];
1565                         value -= duals[iRow];
1566                    }
1567                    for (; j < startPositive_[iSequence+1]; j++) {
1568                         int iRow = indices_[j];
1569                         value += duals[iRow];
1570                    }
1571                    if (value > tolerance) {
1572                         numberWanted--;
1573                         if (value > bestDj) {
1574                              // check flagged variable and correct dj
1575                              if (!model->flagged(iSequence)) {
1576                                   bestDj = value;
1577                                   bestSequence = iSequence;
1578                              } else {
1579                                   // just to make sure we don't exit before got something
1580                                   numberWanted++;
1581                              }
1582                         }
1583                    }
1584                    break;
1585               case ClpSimplex::atLowerBound:
1586                    value = cost[iSequence];
1587                    j = startPositive_[iSequence];
1588                    for (; j < startNegative_[iSequence]; j++) {
1589                         int iRow = indices_[j];
1590                         value -= duals[iRow];
1591                    }
1592                    for (; j < startPositive_[iSequence+1]; j++) {
1593                         int iRow = indices_[j];
1594                         value += duals[iRow];
1595                    }
1596                    value = -value;
1597                    if (value > tolerance) {
1598                         numberWanted--;
1599                         if (value > bestDj) {
1600                              // check flagged variable and correct dj
1601                              if (!model->flagged(iSequence)) {
1602                                   bestDj = value;
1603                                   bestSequence = iSequence;
1604                              } else {
1605                                   // just to make sure we don't exit before got something
1606                                   numberWanted++;
1607                              }
1608                         }
1609                    }
1610                    break;
1611               }
1612          }
1613          if (!numberWanted)
1614               break;
1615     }
1616     if (bestSequence != saveSequence) {
1617          // recompute dj
1618          double value = cost[bestSequence];
1619          j = startPositive_[bestSequence];
1620          for (; j < startNegative_[bestSequence]; j++) {
1621               int iRow = indices_[j];
1622               value -= duals[iRow];
1623          }
1624          for (; j < startPositive_[bestSequence+1]; j++) {
1625               int iRow = indices_[j];
1626               value += duals[iRow];
1627          }
1628          reducedCost[bestSequence] = value;
1629          savedBestSequence_ = bestSequence;
1630          savedBestDj_ = reducedCost[savedBestSequence_];
1631     }
1632     currentWanted_ = numberWanted;
1633}
1634// Allow any parts of a created CoinMatrix to be deleted
1635void
1636ClpPlusMinusOneMatrix::releasePackedMatrix() const
1637{
1638     delete matrix_;
1639     delete [] lengths_;
1640     matrix_ = NULL;
1641     lengths_ = NULL;
1642}
1643/* Returns true if can combine transposeTimes and subsetTransposeTimes
1644   and if it would be faster */
1645bool
1646ClpPlusMinusOneMatrix::canCombine(const ClpSimplex * model,
1647                                  const CoinIndexedVector * pi) const
1648{
1649     int numberInRowArray = pi->getNumElements();
1650     int numberRows = model->numberRows();
1651     bool packed = pi->packedMode();
1652     // factor should be smaller if doing both with two pi vectors
1653     double factor = 0.27;
1654     // We may not want to do by row if there may be cache problems
1655     // It would be nice to find L2 cache size - for moment 512K
1656     // Be slightly optimistic
1657     if (numberColumns_ * sizeof(double) > 1000000) {
1658          if (numberRows * 10 < numberColumns_)
1659               factor *= 0.333333333;
1660          else if (numberRows * 4 < numberColumns_)
1661               factor *= 0.5;
1662          else if (numberRows * 2 < numberColumns_)
1663               factor *= 0.66666666667;
1664          //if (model->numberIterations()%50==0)
1665          //printf("%d nonzero\n",numberInRowArray);
1666     }
1667     // if not packed then bias a bit more towards by column
1668     if (!packed)
1669          factor *= 0.9;
1670     return (numberInRowArray > factor * numberRows || !model->rowCopy());
1671}
1672// These have to match ClpPrimalColumnSteepest version
1673#define reference(i)  (((reference[i>>5]>>(i&31))&1)!=0)
1674// Updates two arrays for steepest
1675void
1676ClpPlusMinusOneMatrix::transposeTimes2(const ClpSimplex * model,
1677                                       const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
1678                                       const CoinIndexedVector * pi2,
1679                                       CoinIndexedVector * spare,
1680                                       double referenceIn, double devex,
1681                                       // Array for exact devex to say what is in reference framework
1682                                       unsigned int * reference,
1683                                       double * weights, double scaleFactor)
1684{
1685     // put row of tableau in dj1
1686     double * pi = pi1->denseVector();
1687     int numberNonZero = 0;
1688     int * index = dj1->getIndices();
1689     double * array = dj1->denseVector();
1690     int numberInRowArray = pi1->getNumElements();
1691     double zeroTolerance = model->zeroTolerance();
1692     bool packed = pi1->packedMode();
1693     // do by column
1694     int iColumn;
1695     assert (!spare->getNumElements());
1696     double * piWeight = pi2->denseVector();
1697     assert (!pi2->packedMode());
1698     bool killDjs = (scaleFactor == 0.0);
1699     if (!scaleFactor)
1700          scaleFactor = 1.0;
1701     // Note scale factor was -1.0
1702     if (packed) {
1703          // need to expand pi into y
1704          assert(spare->capacity() >= model->numberRows());
1705          double * piOld = pi;
1706          pi = spare->denseVector();
1707          const int * whichRow = pi1->getIndices();
1708          int i;
1709          // modify pi so can collapse to one loop
1710          for (i = 0; i < numberInRowArray; i++) {
1711               int iRow = whichRow[i];
1712               pi[iRow] = piOld[i];
1713          }
1714          CoinBigIndex j;
1715          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1716               ClpSimplex::Status status = model->getStatus(iColumn);
1717               if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
1718               double value = 0.0;
1719               for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1720                    int iRow = indices_[j];
1721                    value -= pi[iRow];
1722               }
1723               for (; j < startPositive_[iColumn+1]; j++) {
1724                    int iRow = indices_[j];
1725                    value += pi[iRow];
1726               }
1727               if (fabs(value) > zeroTolerance) {
1728                    // and do other array
1729                    double modification = 0.0;
1730                    for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1731                         int iRow = indices_[j];
1732                         modification += piWeight[iRow];
1733                    }
1734                    for (; j < startPositive_[iColumn+1]; j++) {
1735                         int iRow = indices_[j];
1736                         modification -= piWeight[iRow];
1737                    }
1738                    double thisWeight = weights[iColumn];
1739                    double pivot = value * scaleFactor;
1740                    double pivotSquared = pivot * pivot;
1741                    thisWeight += pivotSquared * devex + pivot * modification;
1742                    if (thisWeight < DEVEX_TRY_NORM) {
1743                         if (referenceIn < 0.0) {
1744                              // steepest
1745                              thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
1746                         } else {
1747                              // exact
1748                              thisWeight = referenceIn * pivotSquared;
1749                              if (reference(iColumn))
1750                                   thisWeight += 1.0;
1751                              thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
1752                         }
1753                    }
1754                    weights[iColumn] = thisWeight;
1755                    if (!killDjs) {
1756                         array[numberNonZero] = value;
1757                         index[numberNonZero++] = iColumn;
1758                    }
1759               }
1760          }
1761          // zero out
1762          for (i = 0; i < numberInRowArray; i++) {
1763               int iRow = whichRow[i];
1764               pi[iRow] = 0.0;
1765          }
1766     } else {
1767          CoinBigIndex j;
1768          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1769               ClpSimplex::Status status = model->getStatus(iColumn);
1770               if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
1771               double value = 0.0;
1772               for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1773                    int iRow = indices_[j];
1774                    value -= pi[iRow];
1775               }
1776               for (; j < startPositive_[iColumn+1]; j++) {
1777                    int iRow = indices_[j];
1778                    value += pi[iRow];
1779               }
1780               if (fabs(value) > zeroTolerance) {
1781                    // and do other array
1782                    double modification = 0.0;
1783                    for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1784                         int iRow = indices_[j];
1785                         modification += piWeight[iRow];
1786                    }
1787                    for (; j < startPositive_[iColumn+1]; j++) {
1788                         int iRow = indices_[j];
1789                         modification -= piWeight[iRow];
1790                    }
1791                    double thisWeight = weights[iColumn];
1792                    double pivot = value * scaleFactor;
1793                    double pivotSquared = pivot * pivot;
1794                    thisWeight += pivotSquared * devex + pivot * modification;
1795                    if (thisWeight < DEVEX_TRY_NORM) {
1796                         if (referenceIn < 0.0) {
1797                              // steepest
1798                              thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
1799                         } else {
1800                              // exact
1801                              thisWeight = referenceIn * pivotSquared;
1802                              if (reference(iColumn))
1803                                   thisWeight += 1.0;
1804                              thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
1805                         }
1806                    }
1807                    weights[iColumn] = thisWeight;
1808                    if (!killDjs) {
1809                         array[iColumn] = value;
1810                         index[numberNonZero++] = iColumn;
1811                    }
1812               }
1813          }
1814     }
1815     dj1->setNumElements(numberNonZero);
1816     spare->setNumElements(0);
1817     if (packed)
1818          dj1->setPackedMode(true);
1819}
1820// Updates second array for steepest and does devex weights
1821void
1822ClpPlusMinusOneMatrix::subsetTimes2(const ClpSimplex * ,
1823                                    CoinIndexedVector * dj1,
1824                                    const CoinIndexedVector * pi2, CoinIndexedVector *,
1825                                    double referenceIn, double devex,
1826                                    // Array for exact devex to say what is in reference framework
1827                                    unsigned int * reference,
1828                                    double * weights, double scaleFactor)
1829{
1830     int number = dj1->getNumElements();
1831     const int * index = dj1->getIndices();
1832     double * array = dj1->denseVector();
1833     assert( dj1->packedMode());
1834
1835     double * piWeight = pi2->denseVector();
1836     bool killDjs = (scaleFactor == 0.0);
1837     if (!scaleFactor)
1838          scaleFactor = 1.0;
1839     for (int k = 0; k < number; k++) {
1840          int iColumn = index[k];
1841          double pivot = array[k] * scaleFactor;
1842          if (killDjs)
1843               array[k] = 0.0;
1844          // and do other array
1845          double modification = 0.0;
1846          CoinBigIndex j;
1847          for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1848               int iRow = indices_[j];
1849               modification += piWeight[iRow];
1850          }
1851          for (; j < startPositive_[iColumn+1]; j++) {
1852               int iRow = indices_[j];
1853               modification -= piWeight[iRow];
1854          }
1855          double thisWeight = weights[iColumn];
1856          double pivotSquared = pivot * pivot;
1857          thisWeight += pivotSquared * devex + pivot * modification;
1858          if (thisWeight < DEVEX_TRY_NORM) {
1859               if (referenceIn < 0.0) {
1860                    // steepest
1861                    thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
1862               } else {
1863                    // exact
1864                    thisWeight = referenceIn * pivotSquared;
1865                    if (reference(iColumn))
1866                         thisWeight += 1.0;
1867                    thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
1868               }
1869          }
1870          weights[iColumn] = thisWeight;
1871     }
1872}
1873/* Set the dimensions of the matrix. In effect, append new empty
1874   columns/rows to the matrix. A negative number for either dimension
1875   means that that dimension doesn't change. Otherwise the new dimensions
1876   MUST be at least as large as the current ones otherwise an exception
1877   is thrown. */
1878void
1879ClpPlusMinusOneMatrix::setDimensions(int newnumrows, int newnumcols)
1880{
1881     if (newnumrows < 0)
1882          newnumrows = numberRows_;
1883     if (newnumrows < numberRows_)
1884          throw CoinError("Bad new rownum (less than current)",
1885                          "setDimensions", "CoinPackedMatrix");
1886
1887     if (newnumcols < 0)
1888          newnumcols = numberColumns_;
1889     if (newnumcols < numberColumns_)
1890          throw CoinError("Bad new colnum (less than current)",
1891                          "setDimensions", "CoinPackedMatrix");
1892
1893     int number = 0;
1894     int length = 0;
1895     if (columnOrdered_) {
1896          length = numberColumns_;
1897          numberColumns_ = newnumcols;
1898          number = numberColumns_;
1899
1900     } else {
1901          length = numberRows_;
1902          numberRows_ = newnumrows;
1903          number = numberRows_;
1904     }
1905     if (number > length) {
1906          CoinBigIndex * temp;
1907          int i;
1908          CoinBigIndex end = startPositive_[length];
1909          temp = new CoinBigIndex [number+1];
1910          CoinMemcpyN(startPositive_, (length + 1), temp);
1911          delete [] startPositive_;
1912          for (i = length + 1; i < number + 1; i++)
1913               temp[i] = end;
1914          startPositive_ = temp;
1915          temp = new CoinBigIndex [number];
1916          CoinMemcpyN(startNegative_, length, temp);
1917          delete [] startNegative_;
1918          for (i = length; i < number; i++)
1919               temp[i] = end;
1920          startNegative_ = temp;
1921     }
1922}
1923#ifndef SLIM_CLP
1924/* Append a set of rows/columns to the end of the matrix. Returns number of errors
1925   i.e. if any of the new rows/columns contain an index that's larger than the
1926   number of columns-1/rows-1 (if numberOther>0) or duplicates
1927   If 0 then rows, 1 if columns */
1928int
1929ClpPlusMinusOneMatrix::appendMatrix(int number, int type,
1930                                    const CoinBigIndex * starts, const int * index,
1931                                    const double * element, int /*numberOther*/)
1932{
1933     int numberErrors = 0;
1934     // make into CoinPackedVector
1935     CoinPackedVectorBase ** vectors =
1936          new CoinPackedVectorBase * [number];
1937     int iVector;
1938     for (iVector = 0; iVector < number; iVector++) {
1939          int iStart = starts[iVector];
1940          vectors[iVector] =
1941               new CoinPackedVector(starts[iVector+1] - iStart,
1942                                    index + iStart, element + iStart);
1943     }
1944     if (type == 0) {
1945          // rows
1946          appendRows(number, vectors);
1947     } else {
1948          // columns
1949          appendCols(number, vectors);
1950     }
1951     for (iVector = 0; iVector < number; iVector++)
1952          delete vectors[iVector];
1953     delete [] vectors;
1954     return numberErrors;
1955}
1956#endif
Note: See TracBrowser for help on using the repository browser.