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

Last change on this file since 800 was 800, checked in by ladanyi, 14 years ago

finishing conversion to svn

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