source: trunk/ClpPlusMinusOneMatrix.cpp @ 284

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

Clp long long stuff

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