source: trunk/ClpPlusMinusOneMatrix.cpp @ 461

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

Trying to make faster

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