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

Last change on this file since 1407 was 1407, checked in by forrest, 12 years ago

fix for ray and warning messages

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