source: trunk/ClpPlusMinusOneMatrix.cpp @ 225

Last change on this file since 225 was 225, checked in by forrest, 16 years ago

This should break everything

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 41.6 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5#include <cstdio>
6
7#include "CoinPragma.hpp"
8#include "CoinIndexedVector.hpp"
9#include "CoinHelperFunctions.hpp"
10
11#include "ClpSimplex.hpp"
12#include "ClpFactorization.hpp"
13// at end to get min/max!
14#include "ClpPlusMinusOneMatrix.hpp"
15#include "ClpMessage.hpp"
16
17//#############################################################################
18// Constructors / Destructor / Assignment
19//#############################################################################
20
21//-------------------------------------------------------------------
22// Default Constructor
23//-------------------------------------------------------------------
24ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix () 
25  : ClpMatrixBase()
26{
27  setType(12);
28  elements_ = NULL;
29  startPositive_ = NULL;
30  startNegative_ = NULL;
31  lengths_=NULL;
32  indices_=NULL;
33  numberRows_=0;
34  numberColumns_=0;
35  columnOrdered_=true;
36}
37
38//-------------------------------------------------------------------
39// Copy constructor
40//-------------------------------------------------------------------
41ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const ClpPlusMinusOneMatrix & rhs) 
42: ClpMatrixBase(rhs)
43{ 
44  elements_ = NULL;
45  startPositive_ = NULL;
46  startNegative_ = NULL;
47  lengths_=NULL;
48  indices_=NULL;
49  numberRows_=rhs.numberRows_;
50  numberColumns_=rhs.numberColumns_;
51  columnOrdered_=rhs.columnOrdered_;
52  if (numberColumns_) {
53    int numberElements = rhs.startPositive_[numberColumns_];
54    indices_ = new int [ numberElements];
55    memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
56    startPositive_ = new int [ numberColumns_+1];
57    memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
58    startNegative_ = new int [ numberColumns_];
59    memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(int));
60  }
61}
62
63ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const CoinPackedMatrix & rhs) 
64  : ClpMatrixBase()
65{ 
66  setType(12);
67  elements_ = NULL;
68  startPositive_ = NULL;
69  startNegative_ = NULL;
70  lengths_=NULL;
71  indices_=NULL;
72  int iColumn;
73  assert (rhs.isColOrdered());
74  // get matrix data pointers
75  const int * row = rhs.getIndices();
76  const CoinBigIndex * columnStart = rhs.getVectorStarts();
77  const int * columnLength = rhs.getVectorLengths(); 
78  const double * elementByColumn = rhs.getElements();
79  numberColumns_ = rhs.getNumCols();
80  bool goodPlusMinusOne=true;
81  numberRows_=-1;
82  indices_ = new int[rhs.getNumElements()];
83  startPositive_ = new int [numberColumns_+1];
84  startNegative_ = new int [numberColumns_];
85  int * temp = new int [rhs.getNumRows()];
86  CoinBigIndex j=0;
87  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
88    CoinBigIndex k;
89    int iNeg=0;
90    startPositive_[iColumn]=j;
91    for (k=columnStart[iColumn];k<columnStart[iColumn]+columnLength[iColumn];
92         k++) {
93      int iRow;
94      if (fabs(elementByColumn[k]-1.0)<1.0e-10) {
95        iRow = row[k];
96        numberRows_ = max(numberRows_,iRow);
97        indices_[j++]=iRow;
98      } else if (fabs(elementByColumn[k]+1.0)<1.0e-10) {
99        iRow = row[k];
100        numberRows_ = max(numberRows_,iRow);
101        temp[iNeg++]=iRow;
102      } else {
103        goodPlusMinusOne = false; // not a network
104      }
105    }
106    if (goodPlusMinusOne) {
107      // move negative
108      startNegative_[iColumn]=j;
109      for (k=0;k<iNeg;k++) {
110        indices_[j++] = temp[k];
111      }
112    } else {
113      break;
114    }
115  }
116  startPositive_[numberColumns_]=j;
117  delete [] temp;
118  if (!goodPlusMinusOne) {
119    delete [] indices_;
120    // put in message
121    indices_=NULL;
122    numberRows_=0;
123    numberColumns_=0;
124    delete [] startPositive_;
125    delete [] startNegative_;
126    startPositive_ = NULL;
127    startNegative_ = NULL;
128  } else {
129    numberRows_ ++; //  correct
130    columnOrdered_ = true;
131  }
132}
133
134//-------------------------------------------------------------------
135// Destructor
136//-------------------------------------------------------------------
137ClpPlusMinusOneMatrix::~ClpPlusMinusOneMatrix ()
138{
139  delete [] elements_;
140  delete [] startPositive_;
141  delete [] startNegative_;
142  delete [] lengths_;
143  delete [] indices_;
144}
145
146//----------------------------------------------------------------
147// Assignment operator
148//-------------------------------------------------------------------
149ClpPlusMinusOneMatrix &
150ClpPlusMinusOneMatrix::operator=(const ClpPlusMinusOneMatrix& rhs)
151{
152  if (this != &rhs) {
153    ClpMatrixBase::operator=(rhs);
154    delete [] elements_;
155    delete [] startPositive_;
156    delete [] startNegative_;
157    delete [] lengths_;
158    delete [] indices_;
159    elements_ = NULL;
160    startPositive_ = NULL;
161    lengths_=NULL;
162    indices_=NULL;
163    numberRows_=rhs.numberRows_;
164    numberColumns_=rhs.numberColumns_;
165    columnOrdered_=rhs.columnOrdered_;
166    if (numberColumns_) {
167      int numberElements = rhs.startPositive_[numberColumns_];
168      indices_ = new int [ numberElements];
169      memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
170      startPositive_ = new int [ numberColumns_+1];
171      memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
172      startNegative_ = new int [ numberColumns_];
173      memcpy(startNegative_,rhs.startNegative_,numberColumns_*sizeof(int));
174    }
175  }
176  return *this;
177}
178//-------------------------------------------------------------------
179// Clone
180//-------------------------------------------------------------------
181ClpMatrixBase * ClpPlusMinusOneMatrix::clone() const
182{
183  return new ClpPlusMinusOneMatrix(*this);
184}
185/* Subset clone (without gaps).  Duplicates are allowed
186   and order is as given */
187ClpMatrixBase * 
188ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
189                              int numberColumns, 
190                              const int * whichColumns) const 
191{
192  return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
193                                   numberColumns, whichColumns);
194}
195/* Subset constructor (without gaps).  Duplicates are allowed
196   and order is as given */
197ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
198                       const ClpPlusMinusOneMatrix & rhs,
199                       int numberRows, const int * whichRow,
200                       int numberColumns, const int * whichColumn)
201: ClpMatrixBase(rhs)
202{ 
203  elements_ = NULL;
204  startPositive_ = NULL;
205  startNegative_ = NULL;
206  lengths_=NULL;
207  indices_=NULL;
208  numberRows_=0;
209  numberColumns_=0;
210  columnOrdered_=rhs.columnOrdered_;
211  if (numberRows<=0||numberColumns<=0) {
212    startPositive_ = new int[1];
213    startPositive_[0] = 0;
214  } else {
215    numberColumns_ = numberColumns;
216    numberRows_ = numberRows;
217    const int * index1 = rhs.indices_;
218    int * startPositive1 = rhs.startPositive_;
219
220    int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
221    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
222    int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
223    int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
224    // Also swap incoming if not column ordered
225    if (!columnOrdered_) {
226      int temp1 = numberRows;
227      numberRows = numberColumns;
228      numberColumns = temp1;
229      const int * temp2;
230      temp2 = whichRow;
231      whichRow = whichColumn;
232      whichColumn = temp2;
233    }
234    // Throw exception if rhs empty
235    if (numberMajor1 <= 0 || numberMinor1 <= 0)
236      throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
237    // Array to say if an old row is in new copy
238    int * newRow = new int [numberMinor1];
239    int iRow;
240    for (iRow=0;iRow<numberMinor1;iRow++) 
241      newRow[iRow] = -1;
242    // and array for duplicating rows
243    int * duplicateRow = new int [numberMinor];
244    int numberBad=0;
245    for (iRow=0;iRow<numberMinor;iRow++) {
246      duplicateRow[iRow] = -1;
247      int kRow = whichRow[iRow];
248      if (kRow>=0  && kRow < numberMinor1) {
249        if (newRow[kRow]<0) {
250          // first time
251          newRow[kRow]=iRow;
252        } else {
253          // duplicate
254          int lastRow = newRow[kRow];
255          newRow[kRow]=iRow;
256          duplicateRow[iRow] = lastRow;
257        }
258      } else {
259        // bad row
260        numberBad++;
261      }
262    }
263
264    if (numberBad)
265      throw CoinError("bad minor entries", 
266                      "subset constructor", "ClpPlusMinusOneMatrix");
267    // now get size and check columns
268    int size = 0;
269    int iColumn;
270    numberBad=0;
271    for (iColumn=0;iColumn<numberMajor;iColumn++) {
272      int kColumn = whichColumn[iColumn];
273      if (kColumn>=0  && kColumn <numberMajor1) {
274        int i;
275        for (i=startPositive1[kColumn];i<startPositive1[kColumn+1];i++) {
276          int kRow = index1[i];
277          kRow = newRow[kRow];
278          while (kRow>=0) {
279            size++;
280            kRow = duplicateRow[kRow];
281          }
282        }
283      } else {
284        // bad column
285        numberBad++;
286        printf("%d %d %d %d\n",iColumn,numberMajor,numberMajor1,kColumn);
287      }
288    }
289    if (numberBad)
290      throw CoinError("bad major entries", 
291                      "subset constructor", "ClpPlusMinusOneMatrix");
292    // now create arrays
293    startPositive_ = new int [numberMajor+1];
294    startNegative_ = new int [numberMajor];
295    indices_ = new int[size];
296    // and fill them
297    size = 0;
298    startPositive_[0]=0;
299    int * startNegative1 = rhs.startNegative_;
300    for (iColumn=0;iColumn<numberMajor;iColumn++) {
301      int kColumn = whichColumn[iColumn];
302      int i;
303      for (i=startPositive1[kColumn];i<startNegative1[kColumn];i++) {
304        int kRow = index1[i];
305        kRow = newRow[kRow];
306        while (kRow>=0) {
307          indices_[size++] = kRow;
308          kRow = duplicateRow[kRow];
309        }
310      }
311      startNegative_[iColumn] = size;
312      for (;i<startPositive1[kColumn+1];i++) {
313        int kRow = index1[i];
314        kRow = newRow[kRow];
315        while (kRow>=0) {
316          indices_[size++] = kRow;
317          kRow = duplicateRow[kRow];
318        }
319      }
320      startPositive_[iColumn+1] = size;
321    }
322    delete [] newRow;
323    delete [] duplicateRow;
324  }
325}
326
327
328/* Returns a new matrix in reverse order without gaps */
329ClpMatrixBase * 
330ClpPlusMinusOneMatrix::reverseOrderedCopy() const
331{
332  int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
333  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
334  // count number in each row/column
335  int * tempP = new int [numberMinor];
336  int * tempN = new int [numberMinor];
337  memset(tempP,0,numberMinor*sizeof(int));
338  memset(tempN,0,numberMinor*sizeof(int));
339  CoinBigIndex j=0;
340  int i;
341  for (i=0;i<numberMajor;i++) {
342    for (;j<startNegative_[i];j++) {
343      int iRow = indices_[j];
344      tempP[iRow]++;
345    }
346    for (;j<startPositive_[i+1];j++) {
347      int iRow = indices_[j];
348      tempN[iRow]++;
349    }
350  }
351  int * newIndices = new int [startPositive_[numberMajor]];
352  int * newP = new int [numberMinor+1];
353  int * newN = new int[numberMinor];
354  int iRow;
355  j=0;
356  // do starts
357  for (iRow=0;iRow<numberMinor;iRow++) {
358    newP[iRow]=j;
359    j += tempP[iRow];
360    tempP[iRow]=newP[iRow];
361    newN[iRow] = j;
362    j += tempN[iRow];
363    tempN[iRow]=newN[iRow];
364  }
365  newP[numberMinor]=j;
366  j=0;
367  for (i=0;i<numberMajor;i++) {
368    for (;j<startNegative_[i];j++) {
369      int iRow = indices_[j];
370      int put = tempP[iRow];
371      newIndices[put++] = i;
372      tempP[iRow] = put;
373    }
374    for (;j<startPositive_[i+1];j++) {
375      int iRow = indices_[j];
376      int put = tempN[iRow];
377      newIndices[put++] = i;
378      tempN[iRow] = put;
379    }
380  }
381  delete [] tempP;
382  delete [] tempN;
383  ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
384  newCopy->passInCopy(numberMinor, numberMajor,
385                      !columnOrdered_,  newIndices, newP, newN);
386  return newCopy;
387}
388//unscaled versions
389void 
390ClpPlusMinusOneMatrix::times(double scalar,
391                   const double * x, double * y) const
392{
393  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
394  int i;
395  CoinBigIndex j;
396  assert (columnOrdered_);
397  for (i=0;i<numberMajor;i++) {
398    double value = scalar*x[i];
399    if (value) {
400      for (j=startPositive_[i];j<startNegative_[i];j++) {
401        int iRow = indices_[j];
402        y[iRow] += value;
403      }
404      for (;j<startPositive_[i+1];j++) {
405        int iRow = indices_[j];
406        y[iRow] -= value;
407      }
408    }
409  }
410}
411void 
412ClpPlusMinusOneMatrix::transposeTimes(double scalar,
413                                const double * x, double * y) const
414{
415  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
416  int i;
417  CoinBigIndex j=0;
418  assert (columnOrdered_);
419  for (i=0;i<numberMajor;i++) {
420    double value = 0.0;
421    for (;j<startNegative_[i];j++) {
422      int iRow = indices_[j];
423      value += x[iRow];
424    }
425    for (;j<startPositive_[i+1];j++) {
426      int iRow = indices_[j];
427      value -= x[iRow];
428    }
429    y[i] += scalar*value;
430  }
431}
432void 
433ClpPlusMinusOneMatrix::times(double scalar,
434                       const double * x, double * y,
435                       const double * rowScale, 
436                       const double * columnScale) const
437{
438  // we know it is not scaled
439  times(scalar, x, y);
440}
441void 
442ClpPlusMinusOneMatrix::transposeTimes( double scalar,
443                                 const double * x, double * y,
444                                 const double * rowScale, 
445                                 const double * columnScale) const
446{
447  // we know it is not scaled
448  transposeTimes(scalar, x, y);
449}
450/* Return <code>x * A + y</code> in <code>z</code>.
451        Squashes small elements and knows about ClpSimplex */
452void 
453ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex * model, double scalar,
454                              const CoinIndexedVector * rowArray,
455                              CoinIndexedVector * y,
456                              CoinIndexedVector * columnArray) const
457{
458  // we know it is not scaled
459  columnArray->clear();
460  double * pi = rowArray->denseVector();
461  int numberNonZero=0;
462  int * index = columnArray->getIndices();
463  double * array = columnArray->denseVector();
464  int numberInRowArray = rowArray->getNumElements();
465  // maybe I need one in OsiSimplex
466  double zeroTolerance = model->factorization()->zeroTolerance();
467  int numberRows = model->numberRows();
468  bool packed = rowArray->packedMode();
469  ClpPlusMinusOneMatrix* rowCopy =
470    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
471  double factor = 0.3;
472  // We may not want to do by row if there may be cache problems
473  int numberColumns = model->numberColumns();
474  // It would be nice to find L2 cache size - for moment 512K
475  // Be slightly optimistic
476  if (numberColumns*sizeof(double)>1000000) {
477    if (numberRows*10<numberColumns)
478      factor=0.1;
479    else if (numberRows*4<numberColumns)
480      factor=0.15;
481    else if (numberRows*2<numberColumns)
482      factor=0.2;
483  }
484  if (numberInRowArray>factor*numberRows||!rowCopy) {
485    assert (!y->getNumElements());
486    // do by column
487    // Need to expand if packed mode
488    int iColumn;
489    CoinBigIndex j=0;
490    assert (columnOrdered_);
491    if (packed) {
492      // need to expand pi into y
493      assert(y->capacity()>=numberRows);
494      double * piOld = pi;
495      pi = y->denseVector();
496      const int * whichRow = rowArray->getIndices();
497      int i;
498      // modify pi so can collapse to one loop
499      for (i=0;i<numberInRowArray;i++) {
500        int iRow = whichRow[i];
501        pi[iRow]=scalar*piOld[i];
502      }
503      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
504        double value = 0.0;
505        for (;j<startNegative_[iColumn];j++) {
506          int iRow = indices_[j];
507          value += pi[iRow];
508        }
509        for (;j<startPositive_[iColumn+1];j++) {
510          int iRow = indices_[j];
511          value -= pi[iRow];
512        }
513        if (fabs(value)>zeroTolerance) {
514          array[numberNonZero]=value;
515          index[numberNonZero++]=iColumn;
516        }
517      }
518      for (i=0;i<numberInRowArray;i++) {
519        int iRow = whichRow[i];
520        pi[iRow]=0.0;
521      }
522    } else {
523      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
524        double value = 0.0;
525        for (;j<startNegative_[iColumn];j++) {
526          int iRow = indices_[j];
527          value += pi[iRow];
528        }
529        for (;j<startPositive_[iColumn+1];j++) {
530          int iRow = indices_[j];
531          value -= pi[iRow];
532        }
533        value *= scalar;
534        if (fabs(value)>zeroTolerance) {
535          index[numberNonZero++]=iColumn;
536          array[iColumn]=value;
537        }
538      }
539    }
540    columnArray->setNumElements(numberNonZero);
541  } else {
542    // do by row
543    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
544  }
545}
546/* Return <code>x * A + y</code> in <code>z</code>.
547        Squashes small elements and knows about ClpSimplex */
548void 
549ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
550                              const CoinIndexedVector * rowArray,
551                              CoinIndexedVector * y,
552                              CoinIndexedVector * columnArray) const
553{
554  columnArray->clear();
555  double * pi = rowArray->denseVector();
556  int numberNonZero=0;
557  int * index = columnArray->getIndices();
558  double * array = columnArray->denseVector();
559  int numberInRowArray = rowArray->getNumElements();
560  // maybe I need one in OsiSimplex
561  double zeroTolerance = model->factorization()->zeroTolerance();
562  const int * column = indices_;
563  const CoinBigIndex * startPositive = startPositive_;
564  const CoinBigIndex * startNegative = startNegative_;
565  const int * whichRow = rowArray->getIndices();
566  bool packed = rowArray->packedMode();
567  if (numberInRowArray>2||y->getNumElements()) {
568    // do by rows
569    int iRow;
570    double * markVector = y->denseVector(); // probably empty .. but
571    int * mark = y->getIndices();
572    int numberOriginal=y->getNumElements();
573    int i;
574    if (packed) {
575      assert(!numberOriginal);
576      numberNonZero=0;
577      // and set up mark as char array
578      char * marked = (char *) (index+columnArray->capacity());
579      double * array2 = y->denseVector();
580#ifdef CLP_DEBUG
581      int numberColumns = model->numberColumns();
582      for (i=0;i<numberColumns;i++) {
583        assert(!marked[i]);
584        assert(!array2[i]);
585      }
586#endif
587      for (i=0;i<numberInRowArray;i++) {
588        iRow = whichRow[i]; 
589        double value = pi[i]*scalar;
590        CoinBigIndex j;
591        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
592          int iColumn = column[j];
593          if (!marked[iColumn]) {
594            marked[iColumn]=1;
595            index[numberNonZero++]=iColumn;
596          }
597          array2[iColumn] += value;
598        }
599        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
600          int iColumn = column[j];
601          if (!marked[iColumn]) {
602            marked[iColumn]=1;
603            index[numberNonZero++]=iColumn;
604          }
605          array2[iColumn] -= value;
606        }
607      }
608      // get rid of tiny values and zero out marked
609      numberOriginal=numberNonZero;
610      numberNonZero=0;
611      for (i=0;i<numberOriginal;i++) {
612        int iColumn = index[i];
613        if (marked[iColumn]) {
614          double value = array2[iColumn];
615          array2[iColumn]=0.0;
616          marked[iColumn]=0;
617          if (fabs(value)>zeroTolerance) {
618            array[numberNonZero]=value;
619            index[numberNonZero++]=iColumn;
620          }
621        }
622      }
623    } else {     
624      for (i=0;i<numberOriginal;i++) {
625        int iColumn = mark[i];
626        index[i]=iColumn;
627        array[iColumn]=markVector[iColumn];
628        markVector[iColumn]=0.0;
629      }
630      numberNonZero=numberOriginal;
631      // and set up mark as char array
632      char * marked = (char *) markVector;
633      for (i=0;i<numberOriginal;i++) {
634        int iColumn = index[i];
635        marked[iColumn]=0;
636      }
637      for (i=0;i<numberInRowArray;i++) {
638        iRow = whichRow[i]; 
639        double value = pi[iRow]*scalar;
640        CoinBigIndex j;
641        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
642          int iColumn = column[j];
643          if (!marked[iColumn]) {
644            marked[iColumn]=1;
645            index[numberNonZero++]=iColumn;
646          }
647          array[iColumn] += value;
648        }
649        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
650          int iColumn = column[j];
651          if (!marked[iColumn]) {
652            marked[iColumn]=1;
653            index[numberNonZero++]=iColumn;
654          }
655          array[iColumn] -= value;
656        }
657      }
658      // get rid of tiny values and zero out marked
659      numberOriginal=numberNonZero;
660      numberNonZero=0;
661      for (i=0;i<numberOriginal;i++) {
662        int iColumn = index[i];
663        marked[iColumn]=0;
664        if (fabs(array[iColumn])>zeroTolerance) {
665          index[numberNonZero++]=iColumn;
666        } else {
667          array[iColumn]=0.0;
668        }
669      }
670    }
671  } else if (numberInRowArray==2) {
672    /* do by rows when two rows (do longer first when not packed
673       and shorter first if packed */
674    int iRow0 = whichRow[0];
675    int iRow1 = whichRow[1];
676    int j;
677    if (packed) {
678      double pi0 = pi[0];
679      double pi1 = pi[1];
680      if (startPositive[iRow0+1]-startPositive[iRow0]>
681          startPositive[iRow1+1]-startPositive[iRow1]) {
682        int temp = iRow0;
683        iRow0 = iRow1;
684        iRow1 = temp;
685        pi0=pi1;
686        pi1=pi[0];
687      }
688      // and set up mark as char array
689      char * marked = (char *) (index+columnArray->capacity());
690      int * lookup = y->getIndices();
691      double value = pi0*scalar;
692      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
693        int iColumn = column[j];
694        array[numberNonZero] = value;
695        marked[iColumn]=1;
696        lookup[iColumn]=numberNonZero;
697        index[numberNonZero++]=iColumn;
698      }
699      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
700        int iColumn = column[j];
701        array[numberNonZero] = -value;
702        marked[iColumn]=1;
703        lookup[iColumn]=numberNonZero;
704        index[numberNonZero++]=iColumn;
705      }
706      int numberOriginal = numberNonZero;
707      value = pi1*scalar;
708      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
709        int iColumn = column[j];
710        if (marked[iColumn]) {
711          int iLookup = lookup[iColumn];
712          array[iLookup] += value;
713        } else {
714          if (fabs(value)>zeroTolerance) {
715            array[numberNonZero] = value;
716            index[numberNonZero++]=iColumn;
717          }
718        }
719      }
720      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
721        int iColumn = column[j];
722        if (marked[iColumn]) {
723          int iLookup = lookup[iColumn];
724          array[iLookup] -= value;
725        } else {
726          if (fabs(value)>zeroTolerance) {
727            array[numberNonZero] = -value;
728            index[numberNonZero++]=iColumn;
729          }
730        }
731      }
732      // get rid of tiny values and zero out marked
733      int nDelete=0;
734      for (j=0;j<numberOriginal;j++) {
735        int iColumn = index[j];
736        marked[iColumn]=0;
737        if (fabs(array[j])<=zeroTolerance) 
738          nDelete++;
739      }
740      if (nDelete) {
741        numberOriginal=numberNonZero;
742        numberNonZero=0;
743        for (j=0;j<numberOriginal;j++) {
744          int iColumn = index[j];
745          double value = array[j];
746          array[j]=0.0;
747          if (fabs(value)>zeroTolerance) {
748            array[numberNonZero]=value;
749            index[numberNonZero++]=iColumn;
750          }
751        }
752      }
753    } else {
754      if (startPositive[iRow0+1]-startPositive[iRow0]<
755          startPositive[iRow1+1]-startPositive[iRow1]) {
756        int temp = iRow0;
757        iRow0 = iRow1;
758        iRow1 = temp;
759      }
760      int numberOriginal;
761      int i;
762      numberNonZero=0;
763      double value;
764      value = pi[iRow0]*scalar;
765      CoinBigIndex j;
766      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
767        int iColumn = column[j];
768        index[numberNonZero++]=iColumn;
769        array[iColumn] = value;
770      }
771      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
772        int iColumn = column[j];
773        index[numberNonZero++]=iColumn;
774        array[iColumn] = -value;
775      }
776      value = pi[iRow1]*scalar;
777      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
778        int iColumn = column[j];
779        double value2= array[iColumn];
780        if (value2) {
781          value2 += value;
782        } else {
783          value2 = value;
784          index[numberNonZero++]=iColumn;
785        }
786        array[iColumn] = value2;
787      }
788      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
789        int iColumn = column[j];
790        double value2= array[iColumn];
791        if (value2) {
792          value2 -= value;
793        } else {
794          value2 = -value;
795          index[numberNonZero++]=iColumn;
796        }
797        array[iColumn] = value2;
798      }
799      // get rid of tiny values and zero out marked
800      numberOriginal=numberNonZero;
801      numberNonZero=0;
802      for (i=0;i<numberOriginal;i++) {
803        int iColumn = index[i];
804        if (fabs(array[iColumn])>zeroTolerance) {
805          index[numberNonZero++]=iColumn;
806        } else {
807          array[iColumn]=0.0;
808        }
809      }
810    }
811  } else if (numberInRowArray==1) {
812    // Just one row
813    int iRow=rowArray->getIndices()[0];
814    numberNonZero=0;
815    double value;
816    iRow = whichRow[0]; 
817    CoinBigIndex j;
818    if (packed) {
819      value = pi[0]*scalar;
820      if (fabs(value)>zeroTolerance) {
821        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
822          int iColumn = column[j];
823          array[numberNonZero] = value;
824          index[numberNonZero++]=iColumn;
825        }
826        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
827          int iColumn = column[j];
828          array[numberNonZero] = -value;
829          index[numberNonZero++]=iColumn;
830        }
831      }
832    } else {
833      value = pi[iRow]*scalar;
834      if (fabs(value)>zeroTolerance) {
835        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
836          int iColumn = column[j];
837          array[iColumn] = value;
838          index[numberNonZero++]=iColumn;
839        }
840        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
841          int iColumn = column[j];
842          array[iColumn] = -value;
843          index[numberNonZero++]=iColumn;
844        }
845      }
846    }
847  }
848  columnArray->setNumElements(numberNonZero);
849  y->setNumElements(0);
850}
851/* Return <code>x *A in <code>z</code> but
852   just for indices in y.
853   Squashes small elements and knows about ClpSimplex */
854void 
855ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex * model,
856                              const CoinIndexedVector * rowArray,
857                              const CoinIndexedVector * y,
858                              CoinIndexedVector * columnArray) const
859{
860  columnArray->clear();
861  double * pi = rowArray->denseVector();
862  int numberNonZero=0;
863  int * index = columnArray->getIndices();
864  double * array = columnArray->denseVector();
865  // maybe I need one in OsiSimplex
866  double zeroTolerance = model->factorization()->zeroTolerance();
867  int jColumn;
868  int numberToDo = y->getNumElements();
869  const int * which = y->getIndices();
870  bool packed = rowArray->packedMode();
871  if (packed) {
872    // need to expand pi into y
873    int numberInRowArray = rowArray->getNumElements();
874    assert(y->capacity()>=model->numberRows());
875    double * piOld = pi;
876    pi = y->denseVector();
877    const int * whichRow = rowArray->getIndices();
878    int i;
879    for (i=0;i<numberInRowArray;i++) {
880      int iRow = whichRow[i];
881      pi[iRow]=piOld[i];
882    }
883    // Must line up with y
884    for (jColumn=0;jColumn<numberToDo;jColumn++) {
885      int iColumn = which[jColumn];
886      double value = 0.0;
887      CoinBigIndex j=startPositive_[iColumn];
888      for (;j<startNegative_[iColumn];j++) {
889        int iRow = indices_[j];
890        value += pi[iRow];
891      }
892      for (;j<startPositive_[iColumn+1];j++) {
893        int iRow = indices_[j];
894        value -= pi[iRow];
895      }
896      array[jColumn]=value;
897    }
898    for (i=0;i<numberInRowArray;i++) {
899      int iRow = whichRow[i];
900      pi[iRow]=0.0;
901    }
902  } else {
903    for (jColumn=0;jColumn<numberToDo;jColumn++) {
904      int iColumn = which[jColumn];
905      double value = 0.0;
906      CoinBigIndex j=startPositive_[iColumn];
907      for (;j<startNegative_[iColumn];j++) {
908        int iRow = indices_[j];
909        value += pi[iRow];
910      }
911      for (;j<startPositive_[iColumn+1];j++) {
912        int iRow = indices_[j];
913        value -= pi[iRow];
914      }
915      if (fabs(value)>zeroTolerance) {
916        index[numberNonZero++]=iColumn;
917        array[iColumn]=value;
918      }
919    }
920  }
921}
922/* Returns number of elements in basis
923   column is basic if entry >=0 */
924CoinBigIndex
925ClpPlusMinusOneMatrix::numberInBasis(const int * columnIsBasic) const 
926{
927  int i;
928  CoinBigIndex numberElements=0;
929  assert (columnOrdered_);
930  for (i=0;i<numberColumns_;i++) {
931    if (columnIsBasic[i]>=0) 
932      numberElements += startPositive_[i+1]-startPositive_[i];
933  }
934  return numberElements;
935}
936// Fills in basis (Returns number of elements and updates numberBasic)
937CoinBigIndex
938ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
939                                const int * columnIsBasic, int & numberBasic,
940                                int * indexRowU, int * indexColumnU,
941                                double * elementU) const 
942{
943#ifdef CLPDEBUG
944  const double * rowScale = model->rowScale();
945  assert (!rowScale);
946#endif
947  int i;
948  CoinBigIndex numberElements=0;
949  assert (columnOrdered_);
950  for (i=0;i<numberColumns_;i++) {
951    if (columnIsBasic[i]>=0) {
952      CoinBigIndex j=startPositive_[i];
953      for (;j<startNegative_[i];j++) {
954        int iRow = indices_[j];
955        indexRowU[numberElements]=iRow;
956        indexColumnU[numberElements]=numberBasic;
957        elementU[numberElements++]=1.0;
958      }
959      for (;j<startPositive_[i+1];j++) {
960        int iRow = indices_[j];
961        indexRowU[numberElements]=iRow;
962        indexColumnU[numberElements]=numberBasic;
963        elementU[numberElements++]=-1.0;
964      }
965      numberBasic++;
966    }
967  }
968  return numberElements;
969}
970/* If element NULL returns number of elements in column part of basis,
971   If not NULL fills in as well */
972CoinBigIndex
973ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
974                                 const int * whichColumn, 
975                                 int numberBasic,
976                                 int numberColumnBasic,
977                                 int * indexRowU, int * indexColumnU,
978                                 double * elementU) const 
979{
980  int i;
981  CoinBigIndex numberElements=0;
982  if (elementU!=NULL) {
983    assert (columnOrdered_);
984    for (i=0;i<numberColumnBasic;i++) {
985      int iColumn = whichColumn[i];
986      CoinBigIndex j=startPositive_[iColumn];
987      for (;j<startNegative_[iColumn];j++) {
988        int iRow = indices_[j];
989        indexRowU[numberElements]=iRow;
990        indexColumnU[numberElements]=numberBasic;
991        elementU[numberElements++]=1.0;
992      }
993      for (;j<startPositive_[iColumn+1];j++) {
994        int iRow = indices_[j];
995        indexRowU[numberElements]=iRow;
996        indexColumnU[numberElements]=numberBasic;
997        elementU[numberElements++]=-1.0;
998      }
999      numberBasic++;
1000    }
1001  } else {
1002    for (i=0;i<numberColumnBasic;i++) {
1003      int iColumn = whichColumn[i];
1004      numberElements += startPositive_[iColumn+1]-startPositive_[iColumn];
1005    }
1006  }
1007  return numberElements;
1008}
1009/* Unpacks a column into an CoinIndexedvector
1010 */
1011void 
1012ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model,
1013                              CoinIndexedVector * rowArray,
1014                              int iColumn) const 
1015{
1016  CoinBigIndex j=startPositive_[iColumn];
1017  for (;j<startNegative_[iColumn];j++) {
1018    int iRow = indices_[j];
1019    rowArray->add(iRow,1.0);
1020  }
1021  for (;j<startPositive_[iColumn+1];j++) {
1022    int iRow = indices_[j];
1023    rowArray->add(iRow,-1.0);
1024  }
1025}
1026/* Unpacks a column into an CoinIndexedvector
1027** in packed foramt
1028Note that model is NOT const.  Bounds and objective could
1029be modified if doing column generation (just for this variable) */
1030void 
1031ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * model,
1032                            CoinIndexedVector * rowArray,
1033                            int iColumn) const
1034{
1035  int * index = rowArray->getIndices();
1036  double * array = rowArray->denseVector();
1037  int number = 0;
1038  CoinBigIndex j=startPositive_[iColumn];
1039  for (;j<startNegative_[iColumn];j++) {
1040    int iRow = indices_[j];
1041    array[number]=1.0;
1042    index[number++]=iRow;
1043  }
1044  for (;j<startPositive_[iColumn+1];j++) {
1045    int iRow = indices_[j];
1046    array[number]=-1.0;
1047    index[number++]=iRow;
1048  }
1049  rowArray->setNumElements(number);
1050  rowArray->setPackedMode(true);
1051}
1052/* Adds multiple of a column into an CoinIndexedvector
1053      You can use quickAdd to add to vector */
1054void 
1055ClpPlusMinusOneMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
1056                   int iColumn, double multiplier) const 
1057{
1058  CoinBigIndex j=startPositive_[iColumn];
1059  for (;j<startNegative_[iColumn];j++) {
1060    int iRow = indices_[j];
1061    rowArray->quickAdd(iRow,multiplier);
1062  }
1063  for (;j<startPositive_[iColumn+1];j++) {
1064    int iRow = indices_[j];
1065    rowArray->quickAdd(iRow,-multiplier);
1066  }
1067}
1068
1069// Return a complete CoinPackedMatrix
1070CoinPackedMatrix * 
1071ClpPlusMinusOneMatrix::getPackedMatrix() const 
1072{
1073  int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
1074  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1075  return new CoinPackedMatrix(columnOrdered_,numberMinor,numberMajor,
1076                              getNumElements(),
1077                              getElements(),indices_,
1078                              startPositive_,getVectorLengths());
1079
1080}
1081/* A vector containing the elements in the packed matrix. Note that there
1082   might be gaps in this list, entries that do not belong to any
1083   major-dimension vector. To get the actual elements one should look at
1084   this vector together with vectorStarts and vectorLengths. */
1085const double * 
1086ClpPlusMinusOneMatrix::getElements() const 
1087{
1088  if (!elements_) {
1089    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1090    int numberElements = startPositive_[numberMajor];
1091    elements_ = new double [numberElements];
1092    CoinBigIndex j=0;
1093    int i;
1094    for (i=0;i<numberMajor;i++) {
1095      for (;j<startNegative_[i];j++) {
1096        elements_[j]=1.0;
1097      }
1098      for (;j<startPositive_[i+1];j++) {
1099        elements_[j]=-1.0;
1100      }
1101    }
1102  }
1103  return elements_;
1104}
1105
1106const CoinBigIndex * 
1107ClpPlusMinusOneMatrix::getVectorStarts() const 
1108{
1109  return startPositive_;
1110}
1111/* The lengths of the major-dimension vectors. */
1112const int * 
1113ClpPlusMinusOneMatrix::getVectorLengths() const
1114{
1115  if (!lengths_) {
1116    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1117    lengths_ = new int [numberMajor];
1118    int i;
1119    for (i=0;i<numberMajor;i++) {
1120      lengths_[i]=startPositive_[i+1]-startPositive_[i];
1121    }
1122  }
1123  return lengths_;
1124}
1125/* Delete the columns whose indices are listed in <code>indDel</code>. */
1126void 
1127ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel) 
1128{
1129  int iColumn;
1130  int newSize=startPositive_[numberColumns_];;
1131  int numberBad=0;
1132  // Use array to make sure we can have duplicates
1133  int * which = new int[numberColumns_];
1134  memset(which,0,numberColumns_*sizeof(int));
1135  int nDuplicate=0;
1136  for (iColumn=0;iColumn<numDel;iColumn++) {
1137    int jColumn = indDel[iColumn];
1138    if (jColumn<0||jColumn>=numberColumns_) {
1139      numberBad++;
1140    } else {
1141      newSize -= startPositive_[jColumn+1]-startPositive_[jColumn];
1142      if (which[jColumn])
1143        nDuplicate++;
1144      else
1145        which[jColumn]=1;
1146    }
1147  }
1148  if (numberBad)
1149    throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
1150  int newNumber = numberColumns_-numDel+nDuplicate;
1151  // Get rid of temporary arrays
1152  delete [] lengths_;
1153  lengths_=NULL;
1154  delete [] elements_;
1155  elements_= NULL;
1156  int * newPositive = new int [newNumber+1];
1157  int * newNegative = new int [newNumber];
1158  int * newIndices = new int [newSize];
1159  newNumber=0;
1160  newSize=0;
1161  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1162    if (!which[iColumn]) {
1163      int start,end;
1164      int i;
1165      start = startPositive_[iColumn];
1166      end=startNegative_[iColumn];
1167      newPositive[newNumber]=newSize;
1168      for (i=start;i<end;i++) 
1169        newIndices[newSize++]=indices_[i];
1170      start = startNegative_[iColumn];
1171      end=startPositive_[iColumn+1];
1172      newNegative[newNumber++]=newSize;
1173      for (i=start;i<end;i++) 
1174        newIndices[newSize++]=indices_[i];
1175    }
1176  }
1177  newPositive[newNumber]=newSize;
1178  delete [] which;
1179  delete [] startPositive_;
1180  startPositive_= newPositive;
1181  delete [] startNegative_;
1182  startNegative_= newNegative;
1183  delete [] indices_;
1184  indices_= newIndices;
1185  numberColumns_ = newNumber;
1186}
1187/* Delete the rows whose indices are listed in <code>indDel</code>. */
1188void 
1189ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel) 
1190{
1191  int iRow;
1192  int numberBad=0;
1193  // Use array to make sure we can have duplicates
1194  int * which = new int[numberRows_];
1195  memset(which,0,numberRows_*sizeof(int));
1196  int nDuplicate=0;
1197  for (iRow=0;iRow<numDel;iRow++) {
1198    int jRow = indDel[iRow];
1199    if (jRow<0||jRow>=numberRows_) {
1200      numberBad++;
1201    } else {
1202      if (which[jRow])
1203        nDuplicate++;
1204      else
1205        which[jRow]=1;
1206    }
1207  }
1208  if (numberBad)
1209    throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
1210  int iElement;
1211  int numberElements=startPositive_[numberColumns_];
1212  int newSize=0;
1213  for (iElement=0;iElement<numberElements;iElement++) {
1214    iRow = indices_[iElement];
1215    if (!which[iRow])
1216      newSize++;
1217  }
1218  int newNumber = numberRows_-numDel+nDuplicate;
1219  // Get rid of temporary arrays
1220  delete [] lengths_;
1221  lengths_=NULL;
1222  delete [] elements_;
1223  elements_= NULL;
1224  int * newIndices = new int [newSize];
1225  newSize=0;
1226  int iColumn;
1227  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1228    int start,end;
1229    int i;
1230    start = startPositive_[iColumn];
1231    end=startNegative_[iColumn];
1232    startPositive_[newNumber]=newSize;
1233    for (i=start;i<end;i++) {
1234        iRow = indices_[i];
1235        if (!which[iRow])
1236          newIndices[newSize++]=iRow;
1237    }
1238    start = startNegative_[iColumn];
1239    end=startPositive_[iColumn+1];
1240    startNegative_[newNumber]=newSize;
1241    for (i=start;i<end;i++) {
1242        iRow = indices_[i];
1243        if (!which[iRow])
1244          newIndices[newSize++]=iRow;
1245    }
1246  }
1247  startPositive_[numberColumns_]=newSize;
1248  delete [] which;
1249  delete [] indices_;
1250  indices_= newIndices;
1251  numberRows_ = newNumber;
1252}
1253bool 
1254ClpPlusMinusOneMatrix::isColOrdered() const 
1255{ 
1256  return columnOrdered_;
1257}
1258/* Number of entries in the packed matrix. */
1259CoinBigIndex
1260ClpPlusMinusOneMatrix::getNumElements() const 
1261{
1262  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1263  if (startPositive_) 
1264    return startPositive_[numberMajor];
1265  else
1266    return 0;
1267}
1268// pass in copy (object takes ownership)
1269void 
1270ClpPlusMinusOneMatrix::passInCopy(int numberRows, int numberColumns,
1271                  bool columnOrdered, int * indices,
1272                  int * startPositive, int * startNegative)
1273{
1274  columnOrdered_=columnOrdered;
1275  startPositive_ = startPositive;
1276  startNegative_ = startNegative;
1277  indices_ = indices;
1278  numberRows_=numberRows;
1279  numberColumns_=numberColumns;
1280}
1281/* Given positive integer weights for each row fills in sum of weights
1282   for each column (and slack).
1283   Returns weights vector
1284*/
1285CoinBigIndex * 
1286ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
1287{
1288  int numberRows = model->numberRows();
1289  int numberColumns =model->numberColumns();
1290  int number = numberRows+numberColumns;
1291  CoinBigIndex * weights = new CoinBigIndex[number];
1292  int i;
1293  for (i=0;i<numberColumns;i++) {
1294    CoinBigIndex j;
1295    CoinBigIndex count=0;
1296    for (j=startPositive_[i];j<startPositive_[i+1];j++) {
1297      int iRow=indices_[j];
1298      count += inputWeights[iRow];
1299    }
1300    weights[i]=count;
1301  }
1302  for (i=0;i<numberRows;i++) {
1303    weights[i+numberColumns]=inputWeights[i];
1304  }
1305  return weights;
1306}
1307// Append Columns
1308void 
1309ClpPlusMinusOneMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
1310{
1311  int iColumn;
1312  int size=0;
1313  int numberBad=0;
1314  for (iColumn=0;iColumn<number;iColumn++) {
1315    int n=columns[iColumn]->getNumElements();
1316    const double * element = columns[iColumn]->getElements();
1317    size += n;
1318    int i;
1319    for (i=0;i<n;i++) {
1320      if (fabs(element[i])!=1.0)
1321        numberBad++;
1322    }
1323  }
1324  if (numberBad)
1325    throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix");
1326  // Get rid of temporary arrays
1327  delete [] lengths_;
1328  lengths_=NULL;
1329  delete [] elements_;
1330  elements_= NULL;
1331  int numberNow = startPositive_[numberColumns_];
1332  int * temp;
1333  temp = new int [numberColumns_+1+number];
1334  memcpy(temp,startPositive_,(numberColumns_+1)*sizeof(int));
1335  delete [] startPositive_;
1336  startPositive_= temp;
1337  temp = new int [numberColumns_+number];
1338  memcpy(temp,startNegative_,numberColumns_*sizeof(int));
1339  delete [] startNegative_;
1340  startNegative_= temp;
1341  temp = new int [numberNow+size];
1342  memcpy(temp,indices_,numberNow*sizeof(int));
1343  delete [] indices_;
1344  indices_= temp;
1345  // now add
1346  size=numberNow;
1347  for (iColumn=0;iColumn<number;iColumn++) {
1348    int n=columns[iColumn]->getNumElements();
1349    const int * row = columns[iColumn]->getIndices();
1350    const double * element = columns[iColumn]->getElements();
1351    int i;
1352    for (i=0;i<n;i++) {
1353      if (element[i]==1.0) 
1354        indices_[size++]=row[i];
1355    }
1356    startNegative_[iColumn+numberColumns_]=size;
1357    for (i=0;i<n;i++) {
1358      if (element[i]==-1.0) 
1359        indices_[size++]=row[i];
1360    }
1361    startPositive_[iColumn+numberColumns_+1]=size;
1362  }
1363 
1364  numberColumns_ += number;
1365}
1366// Append Rows
1367void 
1368ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
1369{
1370  // Allocate arrays to use for counting
1371  int * countPositive = new int [numberColumns_+1];
1372  memset(countPositive,0,numberColumns_*sizeof(int));
1373  int * countNegative = new int [numberColumns_];
1374  memset(countNegative,0,numberColumns_*sizeof(int));
1375  int iRow;
1376  int size=0;
1377  int numberBad=0;
1378  for (iRow=0;iRow<number;iRow++) {
1379    int n=rows[iRow]->getNumElements();
1380    const int * column = rows[iRow]->getIndices();
1381    const double * element = rows[iRow]->getElements();
1382    size += n;
1383    int i;
1384    for (i=0;i<n;i++) {
1385      int iColumn = column[i];
1386      if (element[i]==1.0)
1387        countPositive[iColumn++];
1388      else if (element[i]==-1.0)
1389        countNegative[iColumn++];
1390      else
1391        numberBad++;
1392    }
1393  }
1394  if (numberBad)
1395    throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix");
1396  // Get rid of temporary arrays
1397  delete [] lengths_;
1398  lengths_=NULL;
1399  delete [] elements_;
1400  elements_= NULL;
1401  int numberNow = startPositive_[numberColumns_];
1402  int * newIndices = new int [numberNow+size];
1403  // Update starts and turn counts into positions
1404  // also move current indices
1405  int iColumn;
1406  int numberAdded=0;
1407  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1408    int n,now,move;
1409    now = startPositive_[iColumn];
1410    move = startNegative_[iColumn]-now;
1411    n = countPositive[iColumn];
1412    startPositive_[iColumn] += numberAdded;
1413    memcpy(indices_+now,newIndices+startPositive_[iColumn],move*sizeof(int));
1414    countPositive[iColumn]= startNegative_[iColumn]+numberAdded;
1415    numberAdded += n;
1416    now = startNegative_[iColumn];
1417    move = startPositive_[iColumn+1]-now;
1418    n = countNegative[iColumn];
1419    startNegative_[iColumn] += numberAdded;
1420    memcpy(indices_+now,newIndices+startNegative_[iColumn],move*sizeof(int));
1421    countNegative[iColumn]= startPositive_[iColumn+1]+numberAdded;
1422    numberAdded += n;
1423  }
1424  delete [] indices_;
1425  indices_ = newIndices;
1426  startPositive_[numberColumns_] += numberAdded;
1427  // Now put in
1428  for (iRow=0;iRow<number;iRow++) {
1429    int newRow = numberRows_+iRow;
1430    int n=rows[iRow]->getNumElements();
1431    const int * column = rows[iRow]->getIndices();
1432    const double * element = rows[iRow]->getElements();
1433    int i;
1434    for (i=0;i<n;i++) {
1435      int iColumn = column[i];
1436      int put;
1437      if (element[i]==1.0) {
1438        put = countPositive[iColumn];
1439        countPositive[iColumn] = put+1;
1440      } else {
1441        put = countNegative[iColumn];
1442        countNegative[iColumn] = put+1;
1443      }
1444      indices_[put]=newRow;
1445    }
1446  }
1447  delete [] countPositive;
1448  delete [] countNegative;
1449  numberRows_ += number;
1450}
1451/* Returns largest and smallest elements of both signs.
1452   Largest refers to largest absolute value.
1453*/
1454void 
1455ClpPlusMinusOneMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
1456                       double & smallestPositive, double & largestPositive)
1457{
1458  int iColumn;
1459  bool plusOne=false;
1460  bool minusOne=false;
1461  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
1462    if (startNegative_[iColumn]>startPositive_[iColumn])
1463      plusOne=true;
1464    if (startPositive_[iColumn+1]>startNegative_[iColumn])
1465      minusOne=true;
1466  }
1467  if (minusOne) {
1468    smallestNegative=-1.0;
1469    largestNegative=-1.0;
1470  } else {
1471    smallestNegative=0.0;
1472    largestNegative=0.0;
1473  }
1474  if (plusOne) {
1475    smallestPositive=1.0;
1476    largestPositive=1.0;
1477  } else {
1478    smallestPositive=0.0;
1479    largestPositive=0.0;
1480  }
1481}
Note: See TracBrowser for help on using the repository browser.