source: stable/1.6/Clp/src/ClpPlusMinusOneMatrix.cpp @ 1609

Last change on this file since 1609 was 1111, checked in by forrest, 13 years ago

trying to improve networks

  • 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  matrix_ = 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  matrix_ = 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  matrix_ = 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  matrix_ = 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 matrix_;
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 matrix_;
192    delete [] startPositive_;
193    delete [] startNegative_;
194    delete [] lengths_;
195    delete [] indices_;
196    matrix_ = 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  matrix_ = 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, double * spare) 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  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  memcpy(temp,startPositive_,(numberColumns_+1)*sizeof(CoinBigIndex ));
1346  delete [] startPositive_;
1347  startPositive_= temp;
1348  temp = new CoinBigIndex [numberColumns_+number];
1349  memcpy(temp,startNegative_,numberColumns_*sizeof(CoinBigIndex ));
1350  delete [] startNegative_;
1351  startNegative_= temp;
1352  int * temp2 = new int [numberNow+size];
1353  memcpy(temp2,indices_,numberNow*sizeof(int));
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    memcpy(indices_+now,newIndices+startPositive_[iColumn],move*sizeof(int));
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    memcpy(indices_+now,newIndices+startNegative_[iColumn],move*sizeof(int));
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 = (int) (startFraction*numberColumns_);
1507  int end = CoinMin((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, CoinIndexedVector * dj2,
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->factorization()->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 * model,
1823                              CoinIndexedVector * dj1,
1824                            const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
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    memcpy(temp,startPositive_,(length+1)*sizeof(CoinBigIndex));
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    memcpy(temp,startNegative_,length*sizeof(CoinBigIndex));
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.