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

Last change on this file since 754 was 754, checked in by andreasw, 14 years ago

first version

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