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

Last change on this file since 1879 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

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