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

Last change on this file since 2470 was 2385, checked in by unxusr, 11 months ago

formatting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 110.7 KB
Line 
1/* $Id: ClpPlusMinusOneMatrix.cpp 2385 2019-01-06 19:43:06Z stefan $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6#include <cstdio>
7
8#include "CoinPragma.hpp"
9#include "CoinIndexedVector.hpp"
10#include "CoinPackedVector.hpp"
11#include "CoinHelperFunctions.hpp"
12
13#include "ClpSimplex.hpp"
14#include "ClpFactorization.hpp"
15// at end to get min/max!
16#include "ClpPlusMinusOneMatrix.hpp"
17#include "ClpMessage.hpp"
18#ifdef CLP_PLUS_ONE_MATRIX
19static int oneitcount[13] = { 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0 };
20static void oneit(int i)
21{
22  if (!oneitcount[i])
23    printf("Plus ones for call %d\n", i);
24  oneitcount[i]++;
25  oneitcount[12]++;
26  if ((oneitcount[12] % 1000) == 0) {
27    printf("Plus counts");
28    for (int j = 0; j < 12; j++)
29      printf(" %d", oneitcount[j]);
30    printf("\n");
31  }
32}
33#endif
34//#############################################################################
35// Constructors / Destructor / Assignment
36//#############################################################################
37
38//-------------------------------------------------------------------
39// Default Constructor
40//-------------------------------------------------------------------
41ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix()
42  : ClpMatrixBase()
43{
44  setType(12);
45  matrix_ = NULL;
46  startPositive_ = NULL;
47  startNegative_ = NULL;
48  lengths_ = NULL;
49  indices_ = NULL;
50  numberRows_ = 0;
51  numberColumns_ = 0;
52#ifdef CLP_PLUS_ONE_MATRIX
53  otherFlags_ = 0;
54#endif
55  columnOrdered_ = true;
56}
57
58//-------------------------------------------------------------------
59// Copy constructor
60//-------------------------------------------------------------------
61ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix(const ClpPlusMinusOneMatrix &rhs)
62  : ClpMatrixBase(rhs)
63{
64  matrix_ = NULL;
65  startPositive_ = NULL;
66  startNegative_ = NULL;
67  lengths_ = NULL;
68  indices_ = NULL;
69  numberRows_ = rhs.numberRows_;
70  numberColumns_ = rhs.numberColumns_;
71#ifdef CLP_PLUS_ONE_MATRIX
72  otherFlags_ = rhs.otherFlags_;
73#endif
74  columnOrdered_ = rhs.columnOrdered_;
75  if (numberColumns_) {
76    CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
77    indices_ = new int[numberElements];
78    CoinMemcpyN(rhs.indices_, numberElements, indices_);
79    startPositive_ = new CoinBigIndex[numberColumns_ + 1];
80    CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_);
81    startNegative_ = new CoinBigIndex[numberColumns_];
82    CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_);
83  }
84  int numberRows = getNumRows();
85  if (rhs.rhsOffset_ && numberRows) {
86    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
87  } else {
88    rhsOffset_ = NULL;
89  }
90}
91// Constructor from arrays
92ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix(int numberRows, int numberColumns,
93  bool columnOrdered, const int *indices,
94  const CoinBigIndex *startPositive,
95  const CoinBigIndex *startNegative)
96  : ClpMatrixBase()
97{
98  setType(12);
99  matrix_ = NULL;
100  lengths_ = NULL;
101  numberRows_ = numberRows;
102  numberColumns_ = numberColumns;
103  columnOrdered_ = columnOrdered;
104#ifdef CLP_PLUS_ONE_MATRIX
105  otherFlags_ = 0;
106#endif
107  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
108  CoinBigIndex numberElements = startPositive[numberMajor];
109  startPositive_ = ClpCopyOfArray(startPositive, numberMajor + 1);
110  startNegative_ = ClpCopyOfArray(startNegative, numberMajor);
111  indices_ = ClpCopyOfArray(indices, numberElements);
112  // Check valid
113  checkValid(false);
114}
115
116ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix(const CoinPackedMatrix &rhs)
117  : ClpMatrixBase()
118{
119  setType(12);
120  matrix_ = NULL;
121  startPositive_ = NULL;
122  startNegative_ = NULL;
123  lengths_ = NULL;
124  indices_ = NULL;
125  int iColumn;
126  assert(rhs.isColOrdered());
127  // get matrix data pointers
128  const int *row = rhs.getIndices();
129  const CoinBigIndex *columnStart = rhs.getVectorStarts();
130  const int *columnLength = rhs.getVectorLengths();
131  const double *elementByColumn = rhs.getElements();
132  numberColumns_ = rhs.getNumCols();
133  numberRows_ = -1;
134  indices_ = new int[rhs.getNumElements()];
135  startPositive_ = new CoinBigIndex[numberColumns_ + 1];
136  startNegative_ = new CoinBigIndex[numberColumns_];
137  int *temp = new int[rhs.getNumRows()];
138  CoinBigIndex j = 0;
139  CoinBigIndex numberGoodP = 0;
140  CoinBigIndex numberGoodM = 0;
141  CoinBigIndex numberBad = 0;
142  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
143    CoinBigIndex k;
144    int iNeg = 0;
145    startPositive_[iColumn] = j;
146    for (k = columnStart[iColumn]; k < columnStart[iColumn] + columnLength[iColumn];
147         k++) {
148      int iRow;
149      if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
150        iRow = row[k];
151        numberRows_ = CoinMax(numberRows_, iRow);
152        indices_[j++] = iRow;
153        numberGoodP++;
154      } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
155        iRow = row[k];
156        numberRows_ = CoinMax(numberRows_, iRow);
157        temp[iNeg++] = iRow;
158        numberGoodM++;
159      } else {
160        numberBad++;
161      }
162    }
163    // move negative
164    startNegative_[iColumn] = j;
165    for (k = 0; k < iNeg; k++) {
166      indices_[j++] = temp[k];
167    }
168  }
169  startPositive_[numberColumns_] = j;
170  delete[] temp;
171  if (numberBad) {
172    delete[] indices_;
173    indices_ = NULL;
174    numberRows_ = 0;
175    numberColumns_ = 0;
176    delete[] startPositive_;
177    delete[] startNegative_;
178    // Put in statistics
179    startPositive_ = new CoinBigIndex[3];
180    startPositive_[0] = numberGoodP;
181    startPositive_[1] = numberGoodM;
182    startPositive_[2] = numberBad;
183    startNegative_ = NULL;
184  } else {
185    numberRows_++; //  correct
186    // but number should be same as rhs
187    assert(numberRows_ <= rhs.getNumRows());
188    numberRows_ = rhs.getNumRows();
189    columnOrdered_ = true;
190  }
191  // Check valid
192  if (!numberBad)
193    checkValid(false);
194}
195
196//-------------------------------------------------------------------
197// Destructor
198//-------------------------------------------------------------------
199ClpPlusMinusOneMatrix::~ClpPlusMinusOneMatrix()
200{
201  delete matrix_;
202  delete[] startPositive_;
203  delete[] startNegative_;
204  delete[] lengths_;
205  delete[] indices_;
206}
207
208//----------------------------------------------------------------
209// Assignment operator
210//-------------------------------------------------------------------
211ClpPlusMinusOneMatrix &
212ClpPlusMinusOneMatrix::operator=(const ClpPlusMinusOneMatrix &rhs)
213{
214  if (this != &rhs) {
215    ClpMatrixBase::operator=(rhs);
216    delete matrix_;
217    delete[] startPositive_;
218    delete[] startNegative_;
219    delete[] lengths_;
220    delete[] indices_;
221    matrix_ = NULL;
222    startPositive_ = NULL;
223    lengths_ = NULL;
224    indices_ = NULL;
225    numberRows_ = rhs.numberRows_;
226    numberColumns_ = rhs.numberColumns_;
227    columnOrdered_ = rhs.columnOrdered_;
228#ifdef CLP_PLUS_ONE_MATRIX
229    otherFlags_ = rhs.otherFlags_;
230#endif
231    if (numberColumns_) {
232      CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
233      indices_ = new int[numberElements];
234      CoinMemcpyN(rhs.indices_, numberElements, indices_);
235      startPositive_ = new CoinBigIndex[numberColumns_ + 1];
236      CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_);
237      startNegative_ = new CoinBigIndex[numberColumns_];
238      CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_);
239    }
240  }
241  return *this;
242}
243//-------------------------------------------------------------------
244// Clone
245//-------------------------------------------------------------------
246ClpMatrixBase *ClpPlusMinusOneMatrix::clone() const
247{
248  return new ClpPlusMinusOneMatrix(*this);
249}
250/* Subset clone (without gaps).  Duplicates are allowed
251   and order is as given */
252ClpMatrixBase *
253ClpPlusMinusOneMatrix::subsetClone(int numberRows, const int *whichRows,
254  int numberColumns,
255  const int *whichColumns) const
256{
257  return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
258    numberColumns, whichColumns);
259}
260/* Subset constructor (without gaps).  Duplicates are allowed
261   and order is as given */
262ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix(
263  const ClpPlusMinusOneMatrix &rhs,
264  int numberRows, const int *whichRow,
265  int numberColumns, const int *whichColumn)
266  : ClpMatrixBase(rhs)
267{
268  matrix_ = NULL;
269  startPositive_ = NULL;
270  startNegative_ = NULL;
271  lengths_ = NULL;
272  indices_ = NULL;
273  numberRows_ = 0;
274  numberColumns_ = 0;
275  columnOrdered_ = rhs.columnOrdered_;
276#ifdef CLP_PLUS_ONE_MATRIX
277  otherFlags_ = rhs.otherFlags_;
278#endif
279  if (numberRows <= 0 || numberColumns <= 0) {
280    startPositive_ = new CoinBigIndex[1];
281    startPositive_[0] = 0;
282  } else {
283    numberColumns_ = numberColumns;
284    numberRows_ = numberRows;
285    const int *index1 = rhs.indices_;
286    CoinBigIndex *startPositive1 = rhs.startPositive_;
287
288    int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
289    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
290    int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
291    int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
292    // Also swap incoming if not column ordered
293    if (!columnOrdered_) {
294      int temp1 = numberRows;
295      numberRows = numberColumns;
296      numberColumns = temp1;
297      const int *temp2;
298      temp2 = whichRow;
299      whichRow = whichColumn;
300      whichColumn = temp2;
301    }
302    // Throw exception if rhs empty
303    if (numberMajor1 <= 0 || numberMinor1 <= 0)
304      throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
305    // Array to say if an old row is in new copy
306    int *newRow = new int[numberMinor1];
307    int iRow;
308    for (iRow = 0; iRow < numberMinor1; iRow++)
309      newRow[iRow] = -1;
310    // and array for duplicating rows
311    int *duplicateRow = new int[numberMinor];
312    int numberBad = 0;
313    for (iRow = 0; iRow < numberMinor; iRow++) {
314      duplicateRow[iRow] = -1;
315      int kRow = whichRow[iRow];
316      if (kRow >= 0 && kRow < numberMinor1) {
317        if (newRow[kRow] < 0) {
318          // first time
319          newRow[kRow] = iRow;
320        } else {
321          // duplicate
322          int lastRow = newRow[kRow];
323          newRow[kRow] = iRow;
324          duplicateRow[iRow] = lastRow;
325        }
326      } else {
327        // bad row
328        numberBad++;
329      }
330    }
331
332    if (numberBad)
333      throw CoinError("bad minor entries",
334        "subset constructor", "ClpPlusMinusOneMatrix");
335    // now get size and check columns
336    CoinBigIndex size = 0;
337    int iColumn;
338    numberBad = 0;
339    for (iColumn = 0; iColumn < numberMajor; iColumn++) {
340      int kColumn = whichColumn[iColumn];
341      if (kColumn >= 0 && kColumn < numberMajor1) {
342        CoinBigIndex i;
343        for (i = startPositive1[kColumn]; i < startPositive1[kColumn + 1]; i++) {
344          int kRow = index1[i];
345          kRow = newRow[kRow];
346          while (kRow >= 0) {
347            size++;
348            kRow = duplicateRow[kRow];
349          }
350        }
351      } else {
352        // bad column
353        numberBad++;
354        printf("%d %d %d %d\n", iColumn, numberMajor, numberMajor1, kColumn);
355      }
356    }
357    if (numberBad)
358      throw CoinError("bad major entries",
359        "subset constructor", "ClpPlusMinusOneMatrix");
360    // now create arrays
361    startPositive_ = new CoinBigIndex[numberMajor + 1];
362    startNegative_ = new CoinBigIndex[numberMajor];
363    indices_ = new int[size];
364    // and fill them
365    size = 0;
366    startPositive_[0] = 0;
367    CoinBigIndex *startNegative1 = rhs.startNegative_;
368    for (iColumn = 0; iColumn < numberMajor; iColumn++) {
369      int kColumn = whichColumn[iColumn];
370      CoinBigIndex i;
371      for (i = startPositive1[kColumn]; i < startNegative1[kColumn]; i++) {
372        int kRow = index1[i];
373        kRow = newRow[kRow];
374        while (kRow >= 0) {
375          indices_[size++] = kRow;
376          kRow = duplicateRow[kRow];
377        }
378      }
379      startNegative_[iColumn] = size;
380      for (; i < startPositive1[kColumn + 1]; i++) {
381        int kRow = index1[i];
382        kRow = newRow[kRow];
383        while (kRow >= 0) {
384          indices_[size++] = kRow;
385          kRow = duplicateRow[kRow];
386        }
387      }
388      startPositive_[iColumn + 1] = size;
389    }
390    delete[] newRow;
391    delete[] duplicateRow;
392  }
393  // Check valid
394  checkValid(false);
395}
396
397/* Returns a new matrix in reverse order without gaps */
398ClpMatrixBase *
399ClpPlusMinusOneMatrix::reverseOrderedCopy() const
400{
401  int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
402  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
403  // count number in each row/column
404  CoinBigIndex *tempP = new CoinBigIndex[numberMinor];
405  CoinBigIndex *tempN = new CoinBigIndex[numberMinor];
406  memset(tempP, 0, numberMinor * sizeof(CoinBigIndex));
407  memset(tempN, 0, numberMinor * sizeof(CoinBigIndex));
408  CoinBigIndex j = 0;
409  int i;
410  for (i = 0; i < numberMajor; i++) {
411    for (; j < startNegative_[i]; j++) {
412      int iRow = indices_[j];
413      tempP[iRow]++;
414    }
415    for (; j < startPositive_[i + 1]; j++) {
416      int iRow = indices_[j];
417      tempN[iRow]++;
418    }
419  }
420  int *newIndices = new int[startPositive_[numberMajor]];
421  CoinBigIndex *newP = new CoinBigIndex[numberMinor + 1];
422  CoinBigIndex *newN = new CoinBigIndex[numberMinor];
423  int iRow;
424  j = 0;
425  // do starts
426  for (iRow = 0; iRow < numberMinor; iRow++) {
427    newP[iRow] = j;
428    j += tempP[iRow];
429    tempP[iRow] = newP[iRow];
430    newN[iRow] = j;
431    j += tempN[iRow];
432    tempN[iRow] = newN[iRow];
433  }
434  newP[numberMinor] = j;
435  j = 0;
436  for (i = 0; i < numberMajor; i++) {
437    for (; j < startNegative_[i]; j++) {
438      int iRow = indices_[j];
439      CoinBigIndex put = tempP[iRow];
440      newIndices[put++] = i;
441      tempP[iRow] = put;
442    }
443    for (; j < startPositive_[i + 1]; j++) {
444      int iRow = indices_[j];
445      CoinBigIndex put = tempN[iRow];
446      newIndices[put++] = i;
447      tempN[iRow] = put;
448    }
449  }
450  delete[] tempP;
451  delete[] tempN;
452  ClpPlusMinusOneMatrix *newCopy = new ClpPlusMinusOneMatrix();
453  newCopy->passInCopy(numberMinor, numberMajor,
454    !columnOrdered_, newIndices, newP, newN);
455  return newCopy;
456}
457//static bool doPlusOnes=true;
458//unscaled versions
459void ClpPlusMinusOneMatrix::times(double scalar,
460  const double *COIN_RESTRICT x, double *COIN_RESTRICT y) const
461{
462  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
463  int i;
464  CoinBigIndex j;
465  assert(columnOrdered_);
466#ifdef CLP_PLUS_ONE_MATRIX
467  if ((otherFlags_ & 1) == 0 || !doPlusOnes) {
468#endif
469    for (i = 0; i < numberMajor; i++) {
470      double value = scalar * x[i];
471      if (value) {
472        for (j = startPositive_[i]; j < startNegative_[i]; j++) {
473          int iRow = indices_[j];
474          y[iRow] += value;
475        }
476        for (; j < startPositive_[i + 1]; j++) {
477          int iRow = indices_[j];
478          y[iRow] -= value;
479        }
480      }
481    }
482#ifdef CLP_PLUS_ONE_MATRIX
483  } else {
484    // plus one
485    oneit(0);
486    for (i = 0; i < numberMajor; i++) {
487      double value = scalar * x[i];
488      if (value) {
489        for (j = startPositive_[i]; j < startPositive_[i + 1]; j++) {
490          int iRow = indices_[j];
491          y[iRow] += value;
492        }
493      }
494    }
495  }
496#endif
497}
498void ClpPlusMinusOneMatrix::transposeTimes(double scalar,
499  const double *COIN_RESTRICT x, double *COIN_RESTRICT y) const
500{
501  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
502  int i;
503  CoinBigIndex j = 0;
504  assert(columnOrdered_);
505#ifdef CLP_PLUS_ONE_MATRIX
506  if ((otherFlags_ & 1) == 0 || !doPlusOnes) {
507#endif
508    for (i = 0; i < numberMajor; i++) {
509      double value = 0.0;
510      for (; j < startNegative_[i]; j++) {
511        int iRow = indices_[j];
512        value += x[iRow];
513      }
514      for (; j < startPositive_[i + 1]; j++) {
515        int iRow = indices_[j];
516        value -= x[iRow];
517      }
518      y[i] += scalar * value;
519    }
520#ifdef CLP_PLUS_ONE_MATRIX
521  } else {
522    // plus one
523    oneit(1);
524    for (i = 0; i < numberMajor; i++) {
525      double value = 0.0;
526      for (; j < startPositive_[i + 1]; j++) {
527        int iRow = indices_[j];
528        value += x[iRow];
529      }
530      y[i] += scalar * value;
531    }
532  }
533#endif
534}
535void ClpPlusMinusOneMatrix::times(double scalar,
536  const double *COIN_RESTRICT x, double *COIN_RESTRICT y,
537  const double * /*rowScale*/,
538  const double * /*columnScale*/) const
539{
540  // we know it is not scaled
541  times(scalar, x, y);
542}
543void ClpPlusMinusOneMatrix::transposeTimes(double scalar,
544  const double *COIN_RESTRICT x, double *COIN_RESTRICT y,
545  const double * /*rowScale*/,
546  const double * /*columnScale*/,
547  double * /*spare*/) const
548{
549  // we know it is not scaled
550  transposeTimes(scalar, x, y);
551}
552/* Return <code>x * A + y</code> in <code>z</code>.
553        Squashes small elements and knows about ClpSimplex */
554void ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex *model, double scalar,
555  const CoinIndexedVector *rowArray,
556  CoinIndexedVector *y,
557  CoinIndexedVector *columnArray) const
558{
559  // we know it is not scaled
560  columnArray->clear();
561  double *COIN_RESTRICT pi = rowArray->denseVector();
562  int numberNonZero = 0;
563  int *COIN_RESTRICT index = columnArray->getIndices();
564  double *COIN_RESTRICT array = columnArray->denseVector();
565  int numberInRowArray = rowArray->getNumElements();
566  // maybe I need one in OsiSimplex
567  double zeroTolerance = model->zeroTolerance();
568  int numberRows = model->numberRows();
569  bool packed = rowArray->packedMode();
570#ifndef NO_RTTI
571  ClpPlusMinusOneMatrix *rowCopy = dynamic_cast< ClpPlusMinusOneMatrix * >(model->rowCopy());
572#else
573  ClpPlusMinusOneMatrix *rowCopy = static_cast< ClpPlusMinusOneMatrix * >(model->rowCopy());
574#endif
575  double factor = 0.3;
576  // We may not want to do by row if there may be cache problems
577  int numberColumns = model->numberColumns();
578  // It would be nice to find L2 cache size - for moment 512K
579  // Be slightly optimistic
580  if (numberColumns * sizeof(double) > 1000000) {
581    if (numberRows * 10 < numberColumns)
582      factor = 0.1;
583    else if (numberRows * 4 < numberColumns)
584      factor = 0.15;
585    else if (numberRows * 2 < numberColumns)
586      factor = 0.2;
587  }
588  if (numberInRowArray > factor * numberRows || !rowCopy) {
589    assert(!y->getNumElements());
590    // do by column
591    // Need to expand if packed mode
592    int iColumn;
593    CoinBigIndex j = 0;
594    assert(columnOrdered_);
595    if (packed) {
596      // need to expand pi into y
597      assert(y->capacity() >= numberRows);
598      double *COIN_RESTRICT piOld = pi;
599      pi = y->denseVector();
600      const int *COIN_RESTRICT whichRow = rowArray->getIndices();
601      int i;
602      // modify pi so can collapse to one loop
603      for (i = 0; i < numberInRowArray; i++) {
604        int iRow = whichRow[i];
605        pi[iRow] = scalar * piOld[i];
606      }
607#ifdef CLP_PLUS_ONE_MATRIX
608      if ((otherFlags_ & 1) == 0 || !doPlusOnes) {
609#endif
610        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
611          double value = 0.0;
612          for (; j < startNegative_[iColumn]; j++) {
613            int iRow = indices_[j];
614            value += pi[iRow];
615          }
616          for (; j < startPositive_[iColumn + 1]; j++) {
617            int iRow = indices_[j];
618            value -= pi[iRow];
619          }
620          if (fabs(value) > zeroTolerance) {
621            array[numberNonZero] = value;
622            index[numberNonZero++] = iColumn;
623          }
624        }
625#ifdef CLP_PLUS_ONE_MATRIX
626      } else {
627        // plus one
628        oneit(2);
629        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
630          double value = 0.0;
631          for (; j < startPositive_[iColumn + 1]; j++) {
632            int iRow = indices_[j];
633            value += pi[iRow];
634          }
635          if (fabs(value) > zeroTolerance) {
636            array[numberNonZero] = value;
637            index[numberNonZero++] = iColumn;
638          }
639        }
640      }
641#endif
642      for (i = 0; i < numberInRowArray; i++) {
643        int iRow = whichRow[i];
644        pi[iRow] = 0.0;
645      }
646    } else {
647      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
648        double value = 0.0;
649        for (; j < startNegative_[iColumn]; j++) {
650          int iRow = indices_[j];
651          value += pi[iRow];
652        }
653        for (; j < startPositive_[iColumn + 1]; j++) {
654          int iRow = indices_[j];
655          value -= pi[iRow];
656        }
657        value *= scalar;
658        if (fabs(value) > zeroTolerance) {
659          index[numberNonZero++] = iColumn;
660          array[iColumn] = value;
661        }
662      }
663    }
664    columnArray->setNumElements(numberNonZero);
665  } else {
666    // do by row
667    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
668  }
669}
670/* Return <code>x * A + y</code> in <code>z</code>.
671        Squashes small elements and knows about ClpSimplex */
672void ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar,
673  const CoinIndexedVector *rowArray,
674  CoinIndexedVector *y,
675  CoinIndexedVector *columnArray) const
676{
677  columnArray->clear();
678  double *COIN_RESTRICT pi = rowArray->denseVector();
679  int numberNonZero = 0;
680  int *COIN_RESTRICT index = columnArray->getIndices();
681  double *COIN_RESTRICT array = columnArray->denseVector();
682  int numberInRowArray = rowArray->getNumElements();
683  // maybe I need one in OsiSimplex
684  double zeroTolerance = model->zeroTolerance();
685  const int *COIN_RESTRICT column = indices_;
686  const CoinBigIndex *COIN_RESTRICT startPositive = startPositive_;
687  const CoinBigIndex *COIN_RESTRICT startNegative = startNegative_;
688  const int *COIN_RESTRICT whichRow = rowArray->getIndices();
689  bool packed = rowArray->packedMode();
690  if (numberInRowArray > 2) {
691    // do by rows
692    int iRow;
693    double *COIN_RESTRICT markVector = y->denseVector(); // probably empty .. but
694    int numberOriginal = 0;
695    int i;
696    if (packed) {
697      CoinBigIndex numberCovered = 0;
698      int numberColumns = getNumCols();
699      bool sparse = true;
700      int target = 1 * numberColumns;
701      for (int i = 0; i < numberInRowArray; i++) {
702        int iRow = whichRow[i];
703        numberCovered += startPositive[iRow + 1] - startPositive[iRow];
704        if (numberCovered > target) {
705          sparse = false;
706          break;
707        }
708      }
709      numberNonZero = 0;
710      if (sparse) {
711        // and set up mark as char array
712        char *COIN_RESTRICT marked = reinterpret_cast< char * >(index + columnArray->capacity());
713        double *COIN_RESTRICT array2 = y->denseVector();
714#ifdef CLP_DEBUG
715        int numberColumns = model->numberColumns();
716        for (int i = 0; i < numberColumns; i++) {
717          assert(!marked[i]);
718          assert(!array2[i]);
719        }
720#endif
721#ifdef CLP_PLUS_ONE_MATRIX
722        if ((otherFlags_ & 1) == 0 || !doPlusOnes) {
723#endif
724          for (int i = 0; i < numberInRowArray; i++) {
725            iRow = whichRow[i];
726            double value = pi[i] * scalar;
727            CoinBigIndex j;
728            for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
729              int iColumn = column[j];
730              if (!marked[iColumn]) {
731                marked[iColumn] = 1;
732                index[numberNonZero++] = iColumn;
733              }
734              array2[iColumn] += value;
735            }
736            for (j = startNegative[iRow]; j < startPositive[iRow + 1]; j++) {
737              int iColumn = column[j];
738              if (!marked[iColumn]) {
739                marked[iColumn] = 1;
740                index[numberNonZero++] = iColumn;
741              }
742              array2[iColumn] -= value;
743            }
744          }
745#ifdef CLP_PLUS_ONE_MATRIX
746        } else {
747          // plus one
748          oneit(4);
749          for (int i = 0; i < numberInRowArray; i++) {
750            iRow = whichRow[i];
751            double value = pi[i] * scalar;
752            CoinBigIndex j;
753            for (j = startPositive[iRow]; j < startPositive[iRow + 1]; j++) {
754              int iColumn = column[j];
755              if (!marked[iColumn]) {
756                marked[iColumn] = 1;
757                index[numberNonZero++] = iColumn;
758              }
759              array2[iColumn] += value;
760            }
761          }
762        }
763#endif
764        // get rid of tiny values and zero out marked
765        numberOriginal = numberNonZero;
766        numberNonZero = 0;
767        for (int i = 0; i < numberOriginal; i++) {
768          int iColumn = index[i];
769          if (marked[iColumn]) {
770            double value = array2[iColumn];
771            array2[iColumn] = 0.0;
772            marked[iColumn] = 0;
773            if (fabs(value) > zeroTolerance) {
774              array[numberNonZero] = value;
775              index[numberNonZero++] = iColumn;
776            }
777          }
778        }
779      } else {
780        // not sparse
781#ifdef CLP_PLUS_ONE_MATRIX
782        if ((otherFlags_ & 1) == 0 || !doPlusOnes) {
783#endif
784          for (int i = 0; i < numberInRowArray; i++) {
785            iRow = whichRow[i];
786            double value = pi[i] * scalar;
787            CoinBigIndex j;
788            for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
789              int iColumn = column[j];
790              array[iColumn] += value;
791            }
792            for (j = startNegative[iRow]; j < startPositive[iRow + 1]; j++) {
793              int iColumn = column[j];
794              array[iColumn] -= value;
795            }
796          }
797#ifdef CLP_PLUS_ONE_MATRIX
798        } else {
799          // plus one
800          oneit(5);
801          for (int i = 0; i < numberInRowArray; i++) {
802            iRow = whichRow[i];
803            double value = pi[i] * scalar;
804            CoinBigIndex j;
805            for (j = startPositive[iRow]; j < startPositive[iRow + 1]; j++) {
806              int iColumn = column[j];
807              array[iColumn] += value;
808            }
809          }
810        }
811#endif
812        // get rid of tiny values and count
813        for (int i = 0; i < numberColumns; i++) {
814          double value = array[i];
815          if (value) {
816            array[i] = 0.0;
817            if (fabs(value) > zeroTolerance) {
818              array[numberNonZero] = value;
819              index[numberNonZero++] = i;
820            }
821          }
822        }
823      }
824    } else {
825      numberNonZero = 0;
826      // and set up mark as char array
827      char *COIN_RESTRICT marked = reinterpret_cast< char * >(markVector);
828      for (i = 0; i < numberOriginal; i++) {
829        int iColumn = index[i];
830        marked[iColumn] = 0;
831      }
832      for (i = 0; i < numberInRowArray; i++) {
833        iRow = whichRow[i];
834        double value = pi[iRow] * scalar;
835        CoinBigIndex j;
836        for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
837          int iColumn = column[j];
838          if (!marked[iColumn]) {
839            marked[iColumn] = 1;
840            index[numberNonZero++] = iColumn;
841          }
842          array[iColumn] += value;
843        }
844        for (j = startNegative[iRow]; j < startPositive[iRow + 1]; j++) {
845          int iColumn = column[j];
846          if (!marked[iColumn]) {
847            marked[iColumn] = 1;
848            index[numberNonZero++] = iColumn;
849          }
850          array[iColumn] -= value;
851        }
852      }
853      // get rid of tiny values and zero out marked
854      numberOriginal = numberNonZero;
855      numberNonZero = 0;
856      for (i = 0; i < numberOriginal; i++) {
857        int iColumn = index[i];
858        marked[iColumn] = 0;
859        if (fabs(array[iColumn]) > zeroTolerance) {
860          index[numberNonZero++] = iColumn;
861        } else {
862          array[iColumn] = 0.0;
863        }
864      }
865    }
866  } else if (numberInRowArray == 2) {
867    /* do by rows when two rows (do longer first when not packed
868             and shorter first if packed */
869    int iRow0 = whichRow[0];
870    int iRow1 = whichRow[1];
871    CoinBigIndex j;
872    if (packed) {
873      double pi0 = pi[0];
874      double pi1 = pi[1];
875      if (startPositive[iRow0 + 1] - startPositive[iRow0] > startPositive[iRow1 + 1] - startPositive[iRow1]) {
876        int temp = iRow0;
877        iRow0 = iRow1;
878        iRow1 = temp;
879        pi0 = pi1;
880        pi1 = pi[0];
881      }
882      // and set up mark as char array
883      char *COIN_RESTRICT marked = reinterpret_cast< char * >(index + columnArray->capacity());
884      int *COIN_RESTRICT lookup = y->getIndices();
885      double value = pi0 * scalar;
886      int numberOriginal;
887#ifdef CLP_PLUS_ONE_MATRIX
888      if ((otherFlags_ & 1) == 0 || !doPlusOnes) {
889#endif
890        for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
891          int iColumn = column[j];
892          array[numberNonZero] = value;
893          marked[iColumn] = 1;
894          lookup[iColumn] = numberNonZero;
895          index[numberNonZero++] = iColumn;
896        }
897        for (j = startNegative[iRow0]; j < startPositive[iRow0 + 1]; j++) {
898          int iColumn = column[j];
899          array[numberNonZero] = -value;
900          marked[iColumn] = 1;
901          lookup[iColumn] = numberNonZero;
902          index[numberNonZero++] = iColumn;
903        }
904        numberOriginal = numberNonZero;
905        value = pi1 * scalar;
906        for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
907          int iColumn = column[j];
908          if (marked[iColumn]) {
909            int iLookup = lookup[iColumn];
910            array[iLookup] += value;
911          } else {
912            if (fabs(value) > zeroTolerance) {
913              array[numberNonZero] = value;
914              index[numberNonZero++] = iColumn;
915            }
916          }
917        }
918        for (j = startNegative[iRow1]; j < startPositive[iRow1 + 1]; j++) {
919          int iColumn = column[j];
920          if (marked[iColumn]) {
921            int iLookup = lookup[iColumn];
922            array[iLookup] -= value;
923          } else {
924            if (fabs(value) > zeroTolerance) {
925              array[numberNonZero] = -value;
926              index[numberNonZero++] = iColumn;
927            }
928          }
929        }
930#ifdef CLP_PLUS_ONE_MATRIX
931      } else {
932        // plus one
933        oneit(7);
934        for (j = startPositive[iRow0]; j < startPositive[iRow0 + 1]; j++) {
935          int iColumn = column[j];
936          array[numberNonZero] = value;
937          marked[iColumn] = 1;
938          lookup[iColumn] = numberNonZero;
939          index[numberNonZero++] = iColumn;
940        }
941        numberOriginal = numberNonZero;
942        value = pi1 * scalar;
943        for (j = startPositive[iRow1]; j < startPositive[iRow1 + 1]; j++) {
944          int iColumn = column[j];
945          if (marked[iColumn]) {
946            int iLookup = lookup[iColumn];
947            array[iLookup] += value;
948          } else {
949            if (fabs(value) > zeroTolerance) {
950              array[numberNonZero] = value;
951              index[numberNonZero++] = iColumn;
952            }
953          }
954        }
955      }
956#endif
957      // get rid of tiny values and zero out marked
958      int nDelete = 0;
959      for (j = 0; j < numberOriginal; j++) {
960        int iColumn = index[j];
961        marked[iColumn] = 0;
962        if (fabs(array[j]) <= zeroTolerance)
963          nDelete++;
964      }
965      if (nDelete) {
966        numberOriginal = numberNonZero;
967        numberNonZero = 0;
968        for (j = 0; j < numberOriginal; j++) {
969          int iColumn = index[j];
970          double value = array[j];
971          array[j] = 0.0;
972          if (fabs(value) > zeroTolerance) {
973            array[numberNonZero] = value;
974            index[numberNonZero++] = iColumn;
975          }
976        }
977      }
978    } else {
979      if (startPositive[iRow0 + 1] - startPositive[iRow0] < startPositive[iRow1 + 1] - startPositive[iRow1]) {
980        int temp = iRow0;
981        iRow0 = iRow1;
982        iRow1 = temp;
983      }
984      int numberOriginal;
985      int i;
986      numberNonZero = 0;
987      double value;
988      value = pi[iRow0] * scalar;
989      CoinBigIndex j;
990      for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
991        int iColumn = column[j];
992        index[numberNonZero++] = iColumn;
993        array[iColumn] = value;
994      }
995      for (j = startNegative[iRow0]; j < startPositive[iRow0 + 1]; j++) {
996        int iColumn = column[j];
997        index[numberNonZero++] = iColumn;
998        array[iColumn] = -value;
999      }
1000      value = pi[iRow1] * scalar;
1001      for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
1002        int iColumn = column[j];
1003        double value2 = array[iColumn];
1004        if (value2) {
1005          value2 += value;
1006        } else {
1007          value2 = value;
1008          index[numberNonZero++] = iColumn;
1009        }
1010        array[iColumn] = value2;
1011      }
1012      for (j = startNegative[iRow1]; j < startPositive[iRow1 + 1]; j++) {
1013        int iColumn = column[j];
1014        double value2 = array[iColumn];
1015        if (value2) {
1016          value2 -= value;
1017        } else {
1018          value2 = -value;
1019          index[numberNonZero++] = iColumn;
1020        }
1021        array[iColumn] = value2;
1022      }
1023      // get rid of tiny values and zero out marked
1024      numberOriginal = numberNonZero;
1025      numberNonZero = 0;
1026      for (i = 0; i < numberOriginal; i++) {
1027        int iColumn = index[i];
1028        if (fabs(array[iColumn]) > zeroTolerance) {
1029          index[numberNonZero++] = iColumn;
1030        } else {
1031          array[iColumn] = 0.0;
1032        }
1033      }
1034    }
1035  } else if (numberInRowArray == 1) {
1036    // Just one row
1037    int iRow = rowArray->getIndices()[0];
1038    numberNonZero = 0;
1039    double value;
1040    iRow = whichRow[0];
1041    CoinBigIndex j;
1042    if (packed) {
1043      value = pi[0] * scalar;
1044      if (fabs(value) > zeroTolerance) {
1045#ifdef CLP_PLUS_ONE_MATRIX
1046        if ((otherFlags_ & 1) == 0 || !doPlusOnes) {
1047#endif
1048          for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
1049            int iColumn = column[j];
1050            array[numberNonZero] = value;
1051            index[numberNonZero++] = iColumn;
1052          }
1053          for (j = startNegative[iRow]; j < startPositive[iRow + 1]; j++) {
1054            int iColumn = column[j];
1055            array[numberNonZero] = -value;
1056            index[numberNonZero++] = iColumn;
1057          }
1058#ifdef CLP_PLUS_ONE_MATRIX
1059        } else {
1060          // plus one
1061          oneit(9);
1062          for (j = startPositive[iRow]; j < startPositive[iRow + 1]; j++) {
1063            int iColumn = column[j];
1064            array[numberNonZero] = value;
1065            index[numberNonZero++] = iColumn;
1066          }
1067        }
1068#endif
1069      }
1070    } else {
1071      value = pi[iRow] * scalar;
1072      if (fabs(value) > zeroTolerance) {
1073        for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
1074          int iColumn = column[j];
1075          array[iColumn] = value;
1076          index[numberNonZero++] = iColumn;
1077        }
1078        for (j = startNegative[iRow]; j < startPositive[iRow + 1]; j++) {
1079          int iColumn = column[j];
1080          array[iColumn] = -value;
1081          index[numberNonZero++] = iColumn;
1082        }
1083      }
1084    }
1085  }
1086  columnArray->setNumElements(numberNonZero);
1087  if (packed)
1088    columnArray->setPacked();
1089  y->setNumElements(0);
1090}
1091/* Return <code>x *A in <code>z</code> but
1092   just for indices in y. */
1093void ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex *,
1094  const CoinIndexedVector *rowArray,
1095  const CoinIndexedVector *y,
1096  CoinIndexedVector *columnArray) const
1097{
1098  columnArray->clear();
1099  double *COIN_RESTRICT pi = rowArray->denseVector();
1100  double *COIN_RESTRICT array = columnArray->denseVector();
1101  int jColumn;
1102  int numberToDo = y->getNumElements();
1103  const int *COIN_RESTRICT which = y->getIndices();
1104  assert(!rowArray->packedMode());
1105  columnArray->setPacked();
1106#ifdef CLP_PLUS_ONE_MATRIX
1107  if ((otherFlags_ & 1) == 0 || !doPlusOnes) {
1108#endif
1109    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
1110      int iColumn = which[jColumn];
1111      double value = 0.0;
1112      CoinBigIndex j = startPositive_[iColumn];
1113      for (; j < startNegative_[iColumn]; j++) {
1114        int iRow = indices_[j];
1115        value += pi[iRow];
1116      }
1117      for (; j < startPositive_[iColumn + 1]; j++) {
1118        int iRow = indices_[j];
1119        value -= pi[iRow];
1120      }
1121      array[jColumn] = value;
1122    }
1123#ifdef CLP_PLUS_ONE_MATRIX
1124  } else {
1125    // plus one
1126    oneit(11);
1127    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
1128      int iColumn = which[jColumn];
1129      double value = 0.0;
1130      CoinBigIndex j = startPositive_[iColumn];
1131      for (; j < startPositive_[iColumn + 1]; j++) {
1132        int iRow = indices_[j];
1133        value += pi[iRow];
1134      }
1135      array[jColumn] = value;
1136    }
1137  }
1138#endif
1139}
1140/// returns number of elements in column part of basis,
1141int ClpPlusMinusOneMatrix::countBasis(const int *whichColumn,
1142  int &numberColumnBasic)
1143{
1144  int i;
1145  CoinBigIndex numberElements = 0;
1146  for (i = 0; i < numberColumnBasic; i++) {
1147    int iColumn = whichColumn[i];
1148    numberElements += startPositive_[iColumn + 1] - startPositive_[iColumn];
1149  }
1150#if COIN_BIG_INDEX
1151  if (numberElements > COIN_INT_MAX) {
1152    printf("Factorization too large\n");
1153    abort();
1154  }
1155#endif
1156  return static_cast< int >(numberElements);
1157}
1158void ClpPlusMinusOneMatrix::fillBasis(ClpSimplex *,
1159  const int *whichColumn,
1160  int &numberColumnBasic,
1161  int *indexRowU, int *start,
1162  int *rowCount, int *columnCount,
1163  CoinFactorizationDouble *elementU)
1164{
1165  int i;
1166  CoinBigIndex numberElements = start[0];
1167  assert(columnOrdered_);
1168  for (i = 0; i < numberColumnBasic; i++) {
1169    int iColumn = whichColumn[i];
1170    CoinBigIndex j = startPositive_[iColumn];
1171    for (; j < startNegative_[iColumn]; j++) {
1172      int iRow = indices_[j];
1173      indexRowU[numberElements] = iRow;
1174      rowCount[iRow]++;
1175      elementU[numberElements++] = 1.0;
1176    }
1177    for (; j < startPositive_[iColumn + 1]; j++) {
1178      int iRow = indices_[j];
1179      indexRowU[numberElements] = iRow;
1180      rowCount[iRow]++;
1181      elementU[numberElements++] = -1.0;
1182    }
1183    start[i + 1] = static_cast< int >(numberElements);
1184    columnCount[i] = static_cast< int >(numberElements - start[i]);
1185  }
1186#if COIN_BIG_INDEX
1187  if (numberElements > COIN_INT_MAX)
1188    abort();
1189#endif
1190}
1191/* Unpacks a column into an CoinIndexedvector
1192 */
1193void ClpPlusMinusOneMatrix::unpack(const ClpSimplex *,
1194  CoinIndexedVector *rowArray,
1195  int iColumn) const
1196{
1197  CoinBigIndex j = startPositive_[iColumn];
1198  for (; j < startNegative_[iColumn]; j++) {
1199    int iRow = indices_[j];
1200    rowArray->add(iRow, 1.0);
1201  }
1202  for (; j < startPositive_[iColumn + 1]; j++) {
1203    int iRow = indices_[j];
1204    rowArray->add(iRow, -1.0);
1205  }
1206}
1207/* Unpacks a column into an CoinIndexedvector
1208** in packed foramt
1209Note that model is NOT const.  Bounds and objective could
1210be modified if doing column generation (just for this variable) */
1211void ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex *,
1212  CoinIndexedVector *rowArray,
1213  int iColumn) const
1214{
1215  int *COIN_RESTRICT index = rowArray->getIndices();
1216  double *COIN_RESTRICT array = rowArray->denseVector();
1217  int number = 0;
1218  CoinBigIndex j = startPositive_[iColumn];
1219  for (; j < startNegative_[iColumn]; j++) {
1220    int iRow = indices_[j];
1221    array[number] = 1.0;
1222    index[number++] = iRow;
1223  }
1224  for (; j < startPositive_[iColumn + 1]; j++) {
1225    int iRow = indices_[j];
1226    array[number] = -1.0;
1227    index[number++] = iRow;
1228  }
1229  rowArray->setNumElements(number);
1230  rowArray->setPackedMode(true);
1231}
1232/* Adds multiple of a column into an CoinIndexedvector
1233      You can use quickAdd to add to vector */
1234void ClpPlusMinusOneMatrix::add(const ClpSimplex *, CoinIndexedVector *rowArray,
1235  int iColumn, double multiplier) const
1236{
1237  CoinBigIndex j = startPositive_[iColumn];
1238  for (; j < startNegative_[iColumn]; j++) {
1239    int iRow = indices_[j];
1240    rowArray->quickAdd(iRow, multiplier);
1241  }
1242  for (; j < startPositive_[iColumn + 1]; j++) {
1243    int iRow = indices_[j];
1244    rowArray->quickAdd(iRow, -multiplier);
1245  }
1246}
1247/* Adds multiple of a column into an array */
1248void ClpPlusMinusOneMatrix::add(const ClpSimplex *, double *array,
1249  int iColumn, double multiplier) const
1250{
1251  CoinBigIndex j = startPositive_[iColumn];
1252  for (; j < startNegative_[iColumn]; j++) {
1253    int iRow = indices_[j];
1254    array[iRow] += multiplier;
1255  }
1256  for (; j < startPositive_[iColumn + 1]; j++) {
1257    int iRow = indices_[j];
1258    array[iRow] -= multiplier;
1259  }
1260}
1261
1262// Return a complete CoinPackedMatrix
1263CoinPackedMatrix *
1264ClpPlusMinusOneMatrix::getPackedMatrix() const
1265{
1266  if (!matrix_) {
1267    int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
1268    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1269    CoinBigIndex numberElements = startPositive_[numberMajor];
1270    double *elements = new double[numberElements];
1271    CoinBigIndex j = 0;
1272    int i;
1273    for (i = 0; i < numberMajor; i++) {
1274      for (; j < startNegative_[i]; j++) {
1275        elements[j] = 1.0;
1276      }
1277      for (; j < startPositive_[i + 1]; j++) {
1278        elements[j] = -1.0;
1279      }
1280    }
1281    matrix_ = new CoinPackedMatrix(columnOrdered_, numberMinor, numberMajor,
1282      getNumElements(),
1283      elements, indices_,
1284      startPositive_, getVectorLengths());
1285    delete[] elements;
1286    delete[] lengths_;
1287    lengths_ = NULL;
1288  }
1289  return matrix_;
1290}
1291/* A vector containing the elements in the packed matrix. Note that there
1292   might be gaps in this list, entries that do not belong to any
1293   major-dimension vector. To get the actual elements one should look at
1294   this vector together with vectorStarts and vectorLengths. */
1295const double *
1296ClpPlusMinusOneMatrix::getElements() const
1297{
1298  if (!matrix_)
1299    getPackedMatrix();
1300  return matrix_->getElements();
1301}
1302
1303const CoinBigIndex *
1304ClpPlusMinusOneMatrix::getVectorStarts() const
1305{
1306  return startPositive_;
1307}
1308/* The lengths of the major-dimension vectors. */
1309const int *
1310ClpPlusMinusOneMatrix::getVectorLengths() const
1311{
1312  if (!lengths_) {
1313    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1314    lengths_ = new int[numberMajor];
1315    int i;
1316    for (i = 0; i < numberMajor; i++) {
1317      lengths_[i] = static_cast< int >(startPositive_[i + 1] - startPositive_[i]);
1318    }
1319  }
1320  return lengths_;
1321}
1322/* Delete the columns whose indices are listed in <code>indDel</code>. */
1323void ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int *indDel)
1324{
1325  int iColumn;
1326  CoinBigIndex newSize = startPositive_[numberColumns_];
1327  ;
1328  int numberBad = 0;
1329  // Use array to make sure we can have duplicates
1330  int *which = new int[numberColumns_];
1331  memset(which, 0, numberColumns_ * sizeof(int));
1332  int nDuplicate = 0;
1333  for (iColumn = 0; iColumn < numDel; iColumn++) {
1334    int jColumn = indDel[iColumn];
1335    if (jColumn < 0 || jColumn >= numberColumns_) {
1336      numberBad++;
1337    } else {
1338      newSize -= startPositive_[jColumn + 1] - startPositive_[jColumn];
1339      if (which[jColumn])
1340        nDuplicate++;
1341      else
1342        which[jColumn] = 1;
1343    }
1344  }
1345  if (numberBad)
1346    throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
1347  int newNumber = numberColumns_ - numDel + nDuplicate;
1348  // Get rid of temporary arrays
1349  delete[] lengths_;
1350  lengths_ = NULL;
1351  delete matrix_;
1352  matrix_ = NULL;
1353  CoinBigIndex *newPositive = new CoinBigIndex[newNumber + 1];
1354  CoinBigIndex *newNegative = new CoinBigIndex[newNumber];
1355  int *newIndices = new int[newSize];
1356  newNumber = 0;
1357  newSize = 0;
1358  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1359    if (!which[iColumn]) {
1360      CoinBigIndex start, end;
1361      CoinBigIndex i;
1362      start = startPositive_[iColumn];
1363      end = startNegative_[iColumn];
1364      newPositive[newNumber] = newSize;
1365      for (i = start; i < end; i++)
1366        newIndices[newSize++] = indices_[i];
1367      start = startNegative_[iColumn];
1368      end = startPositive_[iColumn + 1];
1369      newNegative[newNumber++] = newSize;
1370      for (i = start; i < end; i++)
1371        newIndices[newSize++] = indices_[i];
1372    }
1373  }
1374  newPositive[newNumber] = newSize;
1375  delete[] which;
1376  delete[] startPositive_;
1377  startPositive_ = newPositive;
1378  delete[] startNegative_;
1379  startNegative_ = newNegative;
1380  delete[] indices_;
1381  indices_ = newIndices;
1382  numberColumns_ = newNumber;
1383}
1384/* Delete the rows whose indices are listed in <code>indDel</code>. */
1385void ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int *indDel)
1386{
1387  int iRow;
1388  int numberBad = 0;
1389  // Use array to make sure we can have duplicates
1390  int *which = new int[numberRows_];
1391  memset(which, 0, numberRows_ * sizeof(int));
1392  int nDuplicate = 0;
1393  for (iRow = 0; iRow < numDel; iRow++) {
1394    int jRow = indDel[iRow];
1395    if (jRow < 0 || jRow >= numberRows_) {
1396      numberBad++;
1397    } else {
1398      if (which[jRow])
1399        nDuplicate++;
1400      else
1401        which[jRow] = 1;
1402    }
1403  }
1404  if (numberBad)
1405    throw CoinError("Indices out of range", "deleteRows", "ClpPlusMinusOneMatrix");
1406  CoinBigIndex iElement;
1407  CoinBigIndex numberElements = startPositive_[numberColumns_];
1408  CoinBigIndex newSize = 0;
1409  for (iElement = 0; iElement < numberElements; iElement++) {
1410    iRow = indices_[iElement];
1411    if (!which[iRow])
1412      newSize++;
1413  }
1414  int newNumber = numberRows_ - numDel + nDuplicate;
1415  // Get rid of temporary arrays
1416  delete[] lengths_;
1417  lengths_ = NULL;
1418  delete matrix_;
1419  matrix_ = NULL;
1420  // redo which
1421  int numberRows = 0;
1422  for (iRow = 0; iRow < numberRows_; iRow++) {
1423    if (which[iRow]) {
1424      which[iRow] = -1; // not wanted
1425    } else {
1426      which[iRow] = numberRows;
1427      numberRows++;
1428    }
1429  }
1430  int *newIndices = new int[newSize];
1431  newSize = 0;
1432  int iColumn;
1433  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1434    CoinBigIndex start, end;
1435    CoinBigIndex i;
1436    start = startPositive_[iColumn];
1437    end = startNegative_[iColumn];
1438    startPositive_[newNumber] = newSize;
1439    for (i = start; i < end; i++) {
1440      iRow = which[indices_[i]];
1441      if (iRow >= 0)
1442        newIndices[newSize++] = iRow;
1443    }
1444    start = startNegative_[iColumn];
1445    end = startPositive_[iColumn + 1];
1446    startNegative_[newNumber] = newSize;
1447    for (i = start; i < end; i++) {
1448      iRow = which[indices_[i]];
1449      if (iRow >= 0)
1450        newIndices[newSize++] = iRow;
1451    }
1452  }
1453  startPositive_[numberColumns_] = newSize;
1454  delete[] which;
1455  delete[] indices_;
1456  indices_ = newIndices;
1457  numberRows_ = newNumber;
1458}
1459bool ClpPlusMinusOneMatrix::isColOrdered() const
1460{
1461  return columnOrdered_;
1462}
1463/* Number of entries in the packed matrix. */
1464CoinBigIndex
1465ClpPlusMinusOneMatrix::getNumElements() const
1466{
1467  int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1468  if (startPositive_)
1469    return startPositive_[numberMajor];
1470  else
1471    return 0;
1472}
1473// pass in copy (object takes ownership)
1474void ClpPlusMinusOneMatrix::passInCopy(int numberRows, int numberColumns,
1475  bool columnOrdered, int *indices,
1476  CoinBigIndex *startPositive, CoinBigIndex *startNegative)
1477{
1478  columnOrdered_ = columnOrdered;
1479#ifdef CLP_PLUS_ONE_MATRIX
1480  otherFlags_ = 0;
1481#endif
1482  startPositive_ = startPositive;
1483  startNegative_ = startNegative;
1484  indices_ = indices;
1485  numberRows_ = numberRows;
1486  numberColumns_ = numberColumns;
1487  // Check valid
1488  checkValid(false);
1489}
1490// Just checks matrix valid - will say if dimensions not quite right if detail
1491void ClpPlusMinusOneMatrix::checkValid(bool detail) const
1492{
1493  int maxIndex = -1;
1494  int minIndex = columnOrdered_ ? numberRows_ : numberColumns_;
1495  int number = !columnOrdered_ ? numberRows_ : numberColumns_;
1496  CoinBigIndex numberElements = getNumElements();
1497  CoinBigIndex last = -1;
1498  int bad = 0;
1499  // say if all +1
1500#ifdef CLP_PLUS_ONE_MATRIX
1501  otherFlags_ = 1;
1502#endif
1503  for (int i = 0; i < number; i++) {
1504    if (startPositive_[i] < last)
1505      bad++;
1506    else
1507      last = startPositive_[i];
1508#ifdef CLP_PLUS_ONE_MATRIX
1509    if (startNegative_[i] < startPositive_[i + 1])
1510      otherFlags_ &= ~1; // not all +1
1511#endif
1512    if (startNegative_[i] < last)
1513      bad++;
1514    else
1515      last = startNegative_[i];
1516  }
1517  if (startPositive_[number] < last)
1518    bad++;
1519  CoinAssertHint(!bad, "starts are not monotonic");
1520  for (CoinBigIndex cbi = 0; cbi < numberElements; cbi++) {
1521    maxIndex = CoinMax(indices_[cbi], maxIndex);
1522    minIndex = CoinMin(indices_[cbi], minIndex);
1523  }
1524  CoinAssert(maxIndex < (columnOrdered_ ? numberRows_ : numberColumns_));
1525  CoinAssert(minIndex >= 0);
1526  if (detail) {
1527    if (minIndex > 0 || maxIndex + 1 < (columnOrdered_ ? numberRows_ : numberColumns_))
1528      printf("Not full range of indices - %d to %d\n", minIndex, maxIndex);
1529  }
1530}
1531/* Given positive integer weights for each row fills in sum of weights
1532   for each column (and slack).
1533   Returns weights vector
1534*/
1535CoinBigIndex *
1536ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex *model, int *inputWeights) const
1537{
1538  int numberRows = model->numberRows();
1539  int numberColumns = model->numberColumns();
1540  int number = numberRows + numberColumns;
1541  CoinBigIndex *weights = new CoinBigIndex[number];
1542  int i;
1543  for (i = 0; i < numberColumns; i++) {
1544    CoinBigIndex j;
1545    CoinBigIndex count = 0;
1546    for (j = startPositive_[i]; j < startPositive_[i + 1]; j++) {
1547      int iRow = indices_[j];
1548      count += inputWeights[iRow];
1549    }
1550    weights[i] = count;
1551  }
1552  for (i = 0; i < numberRows; i++) {
1553    weights[i + numberColumns] = inputWeights[i];
1554  }
1555  return weights;
1556}
1557// Append Columns
1558void ClpPlusMinusOneMatrix::appendCols(int number, const CoinPackedVectorBase *const *columns)
1559{
1560  int iColumn;
1561  CoinBigIndex size = 0;
1562  int numberBad = 0;
1563  // say if all +1
1564#ifdef CLP_PLUS_ONE_MATRIX
1565  otherFlags_ = 0;
1566#endif
1567  for (iColumn = 0; iColumn < number; iColumn++) {
1568    int n = columns[iColumn]->getNumElements();
1569    const double *element = columns[iColumn]->getElements();
1570    size += n;
1571    int i;
1572    for (i = 0; i < n; i++) {
1573      if (fabs(element[i]) != 1.0)
1574        numberBad++;
1575    }
1576  }
1577  if (numberBad)
1578    throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix");
1579  // Get rid of temporary arrays
1580  delete[] lengths_;
1581  lengths_ = NULL;
1582  delete matrix_;
1583  matrix_ = NULL;
1584  CoinBigIndex numberNow = startPositive_[numberColumns_];
1585  CoinBigIndex *temp;
1586  temp = new CoinBigIndex[numberColumns_ + 1 + number];
1587  CoinMemcpyN(startPositive_, (numberColumns_ + 1), temp);
1588  delete[] startPositive_;
1589  startPositive_ = temp;
1590  temp = new CoinBigIndex[numberColumns_ + number];
1591  CoinMemcpyN(startNegative_, numberColumns_, temp);
1592  delete[] startNegative_;
1593  startNegative_ = temp;
1594  int *temp2 = new int[numberNow + size];
1595  CoinMemcpyN(indices_, numberNow, temp2);
1596  delete[] indices_;
1597  indices_ = temp2;
1598  // now add
1599  size = numberNow;
1600  for (iColumn = 0; iColumn < number; iColumn++) {
1601    int n = columns[iColumn]->getNumElements();
1602    const int *row = columns[iColumn]->getIndices();
1603    const double *element = columns[iColumn]->getElements();
1604    int i;
1605    for (i = 0; i < n; i++) {
1606      if (element[i] == 1.0)
1607        indices_[size++] = row[i];
1608    }
1609    startNegative_[iColumn + numberColumns_] = size;
1610    for (i = 0; i < n; i++) {
1611      if (element[i] == -1.0)
1612        indices_[size++] = row[i];
1613    }
1614    startPositive_[iColumn + numberColumns_ + 1] = size;
1615  }
1616
1617  numberColumns_ += number;
1618}
1619// Append Rows
1620void ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase *const *rows)
1621{
1622  // Allocate arrays to use for counting
1623  int *countPositive = new int[numberColumns_ + 1];
1624  memset(countPositive, 0, numberColumns_ * sizeof(int));
1625  int *countNegative = new int[numberColumns_];
1626  memset(countNegative, 0, numberColumns_ * sizeof(int));
1627  int iRow;
1628  CoinBigIndex size = 0;
1629  int numberBad = 0;
1630  // say if all +1
1631#ifdef CLP_PLUS_ONE_MATRIX
1632  otherFlags_ = 0;
1633#endif
1634  for (iRow = 0; iRow < number; iRow++) {
1635    int n = rows[iRow]->getNumElements();
1636    const int *column = rows[iRow]->getIndices();
1637    const double *element = rows[iRow]->getElements();
1638    size += n;
1639    int i;
1640    for (i = 0; i < n; i++) {
1641      int iColumn = column[i];
1642      if (element[i] == 1.0)
1643        countPositive[iColumn]++;
1644      else if (element[i] == -1.0)
1645        countNegative[iColumn]++;
1646      else
1647        numberBad++;
1648    }
1649  }
1650  if (numberBad)
1651    throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix");
1652  // Get rid of temporary arrays
1653  delete[] lengths_;
1654  lengths_ = NULL;
1655  delete matrix_;
1656  matrix_ = NULL;
1657  CoinBigIndex numberNow = startPositive_[numberColumns_];
1658  int *newIndices = new int[numberNow + size];
1659  // Update starts and turn counts into positions
1660  // also move current indices
1661  int iColumn;
1662  CoinBigIndex numberAdded = 0;
1663  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1664    int n, move;
1665    CoinBigIndex now;
1666    now = startPositive_[iColumn];
1667    move = static_cast< int >(startNegative_[iColumn] - now);
1668    n = countPositive[iColumn];
1669    startPositive_[iColumn] += numberAdded;
1670    CoinMemcpyN(indices_ + now, move, newIndices + startPositive_[iColumn]);
1671    countPositive[iColumn] = static_cast< int >(startNegative_[iColumn] + numberAdded);
1672    numberAdded += n;
1673    now = startNegative_[iColumn];
1674    move = static_cast< int >(startPositive_[iColumn + 1] - now);
1675    n = countNegative[iColumn];
1676    startNegative_[iColumn] += numberAdded;
1677    CoinMemcpyN(indices_ + now, move, newIndices + startNegative_[iColumn]);
1678    countNegative[iColumn] = static_cast< int >(startPositive_[iColumn + 1] + numberAdded);
1679    numberAdded += n;
1680  }
1681  delete[] indices_;
1682  indices_ = newIndices;
1683  startPositive_[numberColumns_] += numberAdded;
1684  // Now put in
1685  for (iRow = 0; iRow < number; iRow++) {
1686    int newRow = numberRows_ + iRow;
1687    int n = rows[iRow]->getNumElements();
1688    const int *column = rows[iRow]->getIndices();
1689    const double *element = rows[iRow]->getElements();
1690    int i;
1691    for (i = 0; i < n; i++) {
1692      int iColumn = column[i];
1693      int put;
1694      if (element[i] == 1.0) {
1695        put = countPositive[iColumn];
1696        countPositive[iColumn] = put + 1;
1697      } else {
1698        put = countNegative[iColumn];
1699        countNegative[iColumn] = put + 1;
1700      }
1701      indices_[put] = newRow;
1702    }
1703  }
1704  delete[] countPositive;
1705  delete[] countNegative;
1706  numberRows_ += number;
1707}
1708/* Returns largest and smallest elements of both signs.
1709   Largest refers to largest absolute value.
1710*/
1711void ClpPlusMinusOneMatrix::rangeOfElements(double &smallestNegative, double &largestNegative,
1712  double &smallestPositive, double &largestPositive)
1713{
1714  int iColumn;
1715  bool plusOne = false;
1716  bool minusOne = false;
1717  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1718    if (startNegative_[iColumn] > startPositive_[iColumn])
1719      plusOne = true;
1720    if (startPositive_[iColumn + 1] > startNegative_[iColumn])
1721      minusOne = true;
1722  }
1723  if (minusOne) {
1724    smallestNegative = -1.0;
1725    largestNegative = -1.0;
1726  } else {
1727    smallestNegative = 0.0;
1728    largestNegative = 0.0;
1729  }
1730  if (plusOne) {
1731    smallestPositive = 1.0;
1732    largestPositive = 1.0;
1733  } else {
1734    smallestPositive = 0.0;
1735    largestPositive = 0.0;
1736  }
1737}
1738// Says whether it can do partial pricing
1739bool ClpPlusMinusOneMatrix::canDoPartialPricing() const
1740{
1741  return true;
1742}
1743// Partial pricing
1744void ClpPlusMinusOneMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
1745  int &bestSequence, int &numberWanted)
1746{
1747  numberWanted = currentWanted_;
1748  int start = static_cast< int >(startFraction * numberColumns_);
1749  int end = CoinMin(static_cast< int >(endFraction * numberColumns_ + 1), numberColumns_);
1750  CoinBigIndex j;
1751  double tolerance = model->currentDualTolerance();
1752  double *COIN_RESTRICT reducedCost = model->djRegion();
1753  const double *COIN_RESTRICT duals = model->dualRowSolution();
1754  const double *COIN_RESTRICT cost = model->costRegion();
1755  double bestDj;
1756  if (bestSequence >= 0)
1757    bestDj = fabs(reducedCost[bestSequence]);
1758  else
1759    bestDj = tolerance;
1760  int sequenceOut = model->sequenceOut();
1761  int saveSequence = bestSequence;
1762  int iSequence;
1763  for (iSequence = start; iSequence < end; iSequence++) {
1764    if (iSequence != sequenceOut) {
1765      double value;
1766      ClpSimplex::Status status = model->getStatus(iSequence);
1767
1768      switch (status) {
1769
1770      case ClpSimplex::basic:
1771      case ClpSimplex::isFixed:
1772        break;
1773      case ClpSimplex::isFree:
1774      case ClpSimplex::superBasic:
1775        value = cost[iSequence];
1776        j = startPositive_[iSequence];
1777        for (; j < startNegative_[iSequence]; j++) {
1778          int iRow = indices_[j];
1779          value -= duals[iRow];
1780        }
1781        for (; j < startPositive_[iSequence + 1]; j++) {
1782          int iRow = indices_[j];
1783          value += duals[iRow];
1784        }
1785        value = fabs(value);
1786        if (value > FREE_ACCEPT * tolerance) {
1787          numberWanted--;
1788          // we are going to bias towards free (but only if reasonable)
1789          value *= FREE_BIAS;
1790          if (value > bestDj) {
1791            // check flagged variable and correct dj
1792            if (!model->flagged(iSequence)) {
1793              bestDj = value;
1794              bestSequence = iSequence;
1795            } else {
1796              // just to make sure we don't exit before got something
1797              numberWanted++;
1798            }
1799          }
1800        }
1801        break;
1802      case ClpSimplex::atUpperBound:
1803        value = cost[iSequence];
1804        j = startPositive_[iSequence];
1805        for (; j < startNegative_[iSequence]; j++) {
1806          int iRow = indices_[j];
1807          value -= duals[iRow];
1808        }
1809        for (; j < startPositive_[iSequence + 1]; j++) {
1810          int iRow = indices_[j];
1811          value += duals[iRow];
1812        }
1813        if (value > tolerance) {
1814          numberWanted--;
1815          if (value > bestDj) {
1816            // check flagged variable and correct dj
1817            if (!model->flagged(iSequence)) {
1818              bestDj = value;
1819              bestSequence = iSequence;
1820            } else {
1821              // just to make sure we don't exit before got something
1822              numberWanted++;
1823            }
1824          }
1825        }
1826        break;
1827      case ClpSimplex::atLowerBound:
1828        value = cost[iSequence];
1829        j = startPositive_[iSequence];
1830        for (; j < startNegative_[iSequence]; j++) {
1831          int iRow = indices_[j];
1832          value -= duals[iRow];
1833        }
1834        for (; j < startPositive_[iSequence + 1]; j++) {
1835          int iRow = indices_[j];
1836          value += duals[iRow];
1837        }
1838        value = -value;
1839        if (value > tolerance) {
1840          numberWanted--;
1841          if (value > bestDj) {
1842            // check flagged variable and correct dj
1843            if (!model->flagged(iSequence)) {
1844              bestDj = value;
1845              bestSequence = iSequence;
1846            } else {
1847              // just to make sure we don't exit before got something
1848              numberWanted++;
1849            }
1850          }
1851        }
1852        break;
1853      }
1854    }
1855    if (!numberWanted)
1856      break;
1857  }
1858  if (bestSequence != saveSequence) {
1859    // recompute dj
1860    double value = cost[bestSequence];
1861    j = startPositive_[bestSequence];
1862    for (; j < startNegative_[bestSequence]; j++) {
1863      int iRow = indices_[j];
1864      value -= duals[iRow];
1865    }
1866    for (; j < startPositive_[bestSequence + 1]; j++) {
1867      int iRow = indices_[j];
1868      value += duals[iRow];
1869    }
1870    reducedCost[bestSequence] = value;
1871    savedBestSequence_ = bestSequence;
1872    savedBestDj_ = reducedCost[savedBestSequence_];
1873  }
1874  currentWanted_ = numberWanted;
1875}
1876// Allow any parts of a created CoinMatrix to be deleted
1877void ClpPlusMinusOneMatrix::releasePackedMatrix() const
1878{
1879  delete matrix_;
1880  delete[] lengths_;
1881  matrix_ = NULL;
1882  lengths_ = NULL;
1883  // say if all +1
1884#ifdef CLP_PLUS_ONE_MATRIX
1885  otherFlags_ = 0;
1886#endif
1887}
1888/* Returns true if can combine transposeTimes and subsetTransposeTimes
1889   and if it would be faster */
1890bool ClpPlusMinusOneMatrix::canCombine(const ClpSimplex *model,
1891  const CoinIndexedVector *pi) const
1892{
1893  int numberInRowArray = pi->getNumElements();
1894  int numberRows = model->numberRows();
1895  bool packed = pi->packedMode();
1896  // factor should be smaller if doing both with two pi vectors
1897  double factor = 0.27;
1898  // We may not want to do by row if there may be cache problems
1899  // It would be nice to find L2 cache size - for moment 512K
1900  // Be slightly optimistic
1901  if (numberColumns_ * sizeof(double) > 1000000) {
1902    if (numberRows * 10 < numberColumns_)
1903      factor *= 0.333333333;
1904    else if (numberRows * 4 < numberColumns_)
1905      factor *= 0.5;
1906    else if (numberRows * 2 < numberColumns_)
1907      factor *= 0.66666666667;
1908    //if (model->numberIterations()%50==0)
1909    //printf("%d nonzero\n",numberInRowArray);
1910  }
1911  // if not packed then bias a bit more towards by column
1912  if (!packed)
1913    factor *= 0.9;
1914  return (numberInRowArray > factor * numberRows || !model->rowCopy());
1915}
1916// These have to match ClpPrimalColumnSteepest version
1917#define reference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0)
1918/* Updates two arrays for steepest and does devex weights
1919   Returns nonzero if updates reduced cost and infeas -
1920   new infeas in dj1  - This does not at present*/
1921int ClpPlusMinusOneMatrix::transposeTimes2(const ClpSimplex *model,
1922  const CoinIndexedVector *pi1, CoinIndexedVector *dj1,
1923  const CoinIndexedVector *pi2,
1924  CoinIndexedVector *spare,
1925  double *, double *,
1926  double referenceIn, double devex,
1927  // Array for exact devex to say what is in reference framework
1928  unsigned int *COIN_RESTRICT reference,
1929  double *COIN_RESTRICT weights, double scaleFactor)
1930{
1931  // put row of tableau in dj1
1932  double *COIN_RESTRICT pi = pi1->denseVector();
1933  int numberNonZero = 0;
1934  int *COIN_RESTRICT index = dj1->getIndices();
1935  double *COIN_RESTRICT array = dj1->denseVector();
1936  int numberInRowArray = pi1->getNumElements();
1937  double zeroTolerance = model->zeroTolerance();
1938  bool packed = pi1->packedMode();
1939  // do by column
1940  int iColumn;
1941  assert(!spare->getNumElements());
1942  double *COIN_RESTRICT piWeight = pi2->denseVector();
1943  assert(!pi2->packedMode());
1944  bool killDjs = (scaleFactor == 0.0);
1945  if (!scaleFactor)
1946    scaleFactor = 1.0;
1947  // Note scale factor was -1.0
1948  if (packed) {
1949    // need to expand pi into y
1950    assert(spare->capacity() >= model->numberRows());
1951    double *COIN_RESTRICT piOld = pi;
1952    pi = spare->denseVector();
1953    const int *COIN_RESTRICT whichRow = pi1->getIndices();
1954    int i;
1955    // modify pi so can collapse to one loop
1956    for (i = 0; i < numberInRowArray; i++) {
1957      int iRow = whichRow[i];
1958      pi[iRow] = piOld[i];
1959    }
1960    CoinBigIndex j;
1961    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1962      ClpSimplex::Status status = model->getStatus(iColumn);
1963      if (status == ClpSimplex::basic || status == ClpSimplex::isFixed)
1964        continue;
1965      double value = 0.0;
1966      for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1967        int iRow = indices_[j];
1968        value -= pi[iRow];
1969      }
1970      for (; j < startPositive_[iColumn + 1]; j++) {
1971        int iRow = indices_[j];
1972        value += pi[iRow];
1973      }
1974      if (fabs(value) > zeroTolerance) {
1975        // and do other array
1976        double modification = 0.0;
1977        for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1978          int iRow = indices_[j];
1979          modification += piWeight[iRow];
1980        }
1981        for (; j < startPositive_[iColumn + 1]; j++) {
1982          int iRow = indices_[j];
1983          modification -= piWeight[iRow];
1984        }
1985        double thisWeight = weights[iColumn];
1986        double pivot = value * scaleFactor;
1987        double pivotSquared = pivot * pivot;
1988        thisWeight += pivotSquared * devex + pivot * modification;
1989        if (thisWeight < DEVEX_TRY_NORM) {
1990          if (referenceIn < 0.0) {
1991            // steepest
1992            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
1993          } else {
1994            // exact
1995            thisWeight = referenceIn * pivotSquared;
1996            if (reference(iColumn))
1997              thisWeight += 1.0;
1998            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
1999          }
2000        }
2001        weights[iColumn] = thisWeight;
2002        if (!killDjs) {
2003          array[numberNonZero] = value;
2004          index[numberNonZero++] = iColumn;
2005        }
2006      }
2007    }
2008    // zero out
2009    for (i = 0; i < numberInRowArray; i++) {
2010      int iRow = whichRow[i];
2011      pi[iRow] = 0.0;
2012    }
2013  } else {
2014    CoinBigIndex j;
2015    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2016      ClpSimplex::Status status = model->getStatus(iColumn);
2017      if (status == ClpSimplex::basic || status == ClpSimplex::isFixed)
2018        continue;
2019      double value = 0.0;
2020      for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
2021        int iRow = indices_[j];
2022        value -= pi[iRow];
2023      }
2024      for (; j < startPositive_[iColumn + 1]; j++) {
2025        int iRow = indices_[j];
2026        value += pi[iRow];
2027      }
2028      if (fabs(value) > zeroTolerance) {
2029        // and do other array
2030        double modification = 0.0;
2031        for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
2032          int iRow = indices_[j];
2033          modification += piWeight[iRow];
2034        }
2035        for (; j < startPositive_[iColumn + 1]; j++) {
2036          int iRow = indices_[j];
2037          modification -= piWeight[iRow];
2038        }
2039        double thisWeight = weights[iColumn];
2040        double pivot = value * scaleFactor;
2041        double pivotSquared = pivot * pivot;
2042        thisWeight += pivotSquared * devex + pivot * modification;
2043        if (thisWeight < DEVEX_TRY_NORM) {
2044          if (referenceIn < 0.0) {
2045            // steepest
2046            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2047          } else {
2048            // exact
2049            thisWeight = referenceIn * pivotSquared;
2050            if (reference(iColumn))
2051              thisWeight += 1.0;
2052            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2053          }
2054        }
2055        weights[iColumn] = thisWeight;
2056        if (!killDjs) {
2057          array[iColumn] = value;
2058          index[numberNonZero++] = iColumn;
2059        }
2060      }
2061    }
2062  }
2063  dj1->setNumElements(numberNonZero);
2064  spare->setNumElements(0);
2065  if (packed)
2066    dj1->setPackedMode(true);
2067  return 0;
2068}
2069// Updates second array for steepest and does devex weights
2070void ClpPlusMinusOneMatrix::subsetTimes2(const ClpSimplex *,
2071  CoinIndexedVector *dj1,
2072  const CoinIndexedVector *pi2, CoinIndexedVector *,
2073  double referenceIn, double devex,
2074  // Array for exact devex to say what is in reference framework
2075  unsigned int *COIN_RESTRICT reference,
2076  double *COIN_RESTRICT weights, double scaleFactor)
2077{
2078  int number = dj1->getNumElements();
2079  const int *COIN_RESTRICT index = dj1->getIndices();
2080  double *COIN_RESTRICT array = dj1->denseVector();
2081  assert(dj1->packedMode());
2082
2083  double *COIN_RESTRICT piWeight = pi2->denseVector();
2084  bool killDjs = (scaleFactor == 0.0);
2085  if (!scaleFactor)
2086    scaleFactor = 1.0;
2087  for (int k = 0; k < number; k++) {
2088    int iColumn = index[k];
2089    double pivot = array[k] * scaleFactor;
2090    if (killDjs)
2091      array[k] = 0.0;
2092    // and do other array
2093    double modification = 0.0;
2094    CoinBigIndex j;
2095    for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
2096      int iRow = indices_[j];
2097      modification += piWeight[iRow];
2098    }
2099    for (; j < startPositive_[iColumn + 1]; j++) {
2100      int iRow = indices_[j];
2101      modification -= piWeight[iRow];
2102    }
2103    double thisWeight = weights[iColumn];
2104    double pivotSquared = pivot * pivot;
2105    thisWeight += pivotSquared * devex + pivot * modification;
2106    if (thisWeight < DEVEX_TRY_NORM) {
2107      if (referenceIn < 0.0) {
2108        // steepest
2109        thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2110      } else {
2111        // exact
2112        thisWeight = referenceIn * pivotSquared;
2113        if (reference(iColumn))
2114          thisWeight += 1.0;
2115        thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2116      }
2117    }
2118    weights[iColumn] = thisWeight;
2119  }
2120}
2121/* Set the dimensions of the matrix. In effect, append new empty
2122   columns/rows to the matrix. A negative number for either dimension
2123   means that that dimension doesn't change. Otherwise the new dimensions
2124   MUST be at least as large as the current ones otherwise an exception
2125   is thrown. */
2126void ClpPlusMinusOneMatrix::setDimensions(int newnumrows, int newnumcols)
2127{
2128  if (newnumrows < 0)
2129    newnumrows = numberRows_;
2130  if (newnumrows < numberRows_)
2131    throw CoinError("Bad new rownum (less than current)",
2132      "setDimensions", "CoinPackedMatrix");
2133
2134  if (newnumcols < 0)
2135    newnumcols = numberColumns_;
2136  if (newnumcols < numberColumns_)
2137    throw CoinError("Bad new colnum (less than current)",
2138      "setDimensions", "CoinPackedMatrix");
2139
2140  int number = 0;
2141  int length = 0;
2142  if (columnOrdered_) {
2143    length = numberColumns_;
2144    numberColumns_ = newnumcols;
2145    number = numberColumns_;
2146
2147  } else {
2148    length = numberRows_;
2149    numberRows_ = newnumrows;
2150    number = numberRows_;
2151  }
2152  if (number > length) {
2153    CoinBigIndex *temp;
2154    int i;
2155    CoinBigIndex end = startPositive_[length];
2156    temp = new CoinBigIndex[number + 1];
2157    CoinMemcpyN(startPositive_, (length + 1), temp);
2158    delete[] startPositive_;
2159    for (i = length + 1; i < number + 1; i++)
2160      temp[i] = end;
2161    startPositive_ = temp;
2162    temp = new CoinBigIndex[number];
2163    CoinMemcpyN(startNegative_, length, temp);
2164    delete[] startNegative_;
2165    for (i = length; i < number; i++)
2166      temp[i] = end;
2167    startNegative_ = temp;
2168  }
2169}
2170#ifndef SLIM_CLP
2171/* Append a set of rows/columns to the end of the matrix. Returns number of errors
2172   i.e. if any of the new rows/columns contain an index that's larger than the
2173   number of columns-1/rows-1 (if numberOther>0) or duplicates
2174   If 0 then rows, 1 if columns */
2175int ClpPlusMinusOneMatrix::appendMatrix(int number, int type,
2176  const CoinBigIndex *starts, const int *index,
2177  const double *element, int /*numberOther*/)
2178{
2179  int numberErrors = 0;
2180  // say if all +1
2181#ifdef CLP_PLUS_ONE_MATRIX
2182  otherFlags_ = 0;
2183#endif
2184  // make into CoinPackedVector
2185  CoinPackedVectorBase **vectors = new CoinPackedVectorBase *[number];
2186  int iVector;
2187  for (iVector = 0; iVector < number; iVector++) {
2188    CoinBigIndex iStart = starts[iVector];
2189    vectors[iVector] = new CoinPackedVector(static_cast< int >(starts[iVector + 1] - iStart),
2190      index + iStart, element + iStart);
2191  }
2192  if (type == 0) {
2193    // rows
2194    appendRows(number, vectors);
2195  } else {
2196    // columns
2197    appendCols(number, vectors);
2198  }
2199  for (iVector = 0; iVector < number; iVector++)
2200    delete vectors[iVector];
2201  delete[] vectors;
2202  return numberErrors;
2203}
2204#endif
2205#if CLP_POOL_MATRIX
2206//#############################################################################
2207// Constructors / Destructor / Assignment
2208//#############################################################################
2209
2210//-------------------------------------------------------------------
2211// Default Constructor
2212//-------------------------------------------------------------------
2213ClpPoolMatrix::ClpPoolMatrix()
2214  : ClpMatrixBase()
2215{
2216  setType(13);
2217  matrix_ = NULL;
2218  lengths_ = NULL;
2219  elements_ = NULL;
2220  columnStart_ = NULL;
2221  stuff_ = NULL;
2222  numberRows_ = 0;
2223  numberColumns_ = 0;
2224  numberDifferent_ = 0;
2225}
2226
2227//-------------------------------------------------------------------
2228// Copy constructor
2229//-------------------------------------------------------------------
2230ClpPoolMatrix::ClpPoolMatrix(const ClpPoolMatrix &rhs)
2231  : ClpMatrixBase(rhs)
2232{
2233  setType(13);
2234  matrix_ = NULL;
2235  lengths_ = NULL;
2236  elements_ = NULL;
2237  columnStart_ = NULL;
2238  stuff_ = NULL;
2239  numberRows_ = rhs.numberRows_;
2240  numberColumns_ = rhs.numberColumns_;
2241  numberDifferent_ = rhs.numberDifferent_;
2242  if (numberColumns_) {
2243    columnStart_ = CoinCopyOfArray(rhs.columnStart_, numberColumns_ + 1);
2244    CoinBigIndex numberElements = columnStart_[numberColumns_];
2245    stuff_ = CoinCopyOfArray(rhs.stuff_, numberElements);
2246    elements_ = CoinCopyOfArray(rhs.elements_, numberDifferent_);
2247  }
2248}
2249static const int mmult[] = {
2250  262139, 259459, 256889, 254291, 251701, 249133, 246709, 244247
2251};
2252inline int hashit(double value)
2253{
2254  const unsigned char *chars = reinterpret_cast< unsigned char * >(&value);
2255  int n = 0;
2256  for (int j = 0; j < 8; j++)
2257    n += mmult[j] * chars[j];
2258  n = n % (1 << CLP_POOL_SIZE);
2259  return n;
2260}
2261ClpPoolMatrix::ClpPoolMatrix(const CoinPackedMatrix &rhs)
2262  : ClpMatrixBase()
2263{
2264  setType(13);
2265  matrix_ = NULL;
2266  lengths_ = NULL;
2267  numberDifferent_ = 0;
2268  assert(rhs.isColOrdered());
2269  // get matrix data pointers
2270  const int *row = rhs.getIndices();
2271  const CoinBigIndex *columnStart = rhs.getVectorStarts();
2272  const int *columnLength = rhs.getVectorLengths();
2273  const double *elementByColumn = rhs.getElements();
2274  numberColumns_ = rhs.getNumCols();
2275  numberRows_ = rhs.getNumRows();
2276  assert(numberRows_ < (1 << CLP_POOL_MATRIX));
2277  elements_ = NULL;
2278  columnStart_ = new CoinBigIndex[numberColumns_ + 1];
2279  stuff_ = new poolInfo[rhs.getNumElements()];
2280  int maxPool = 1 << CLP_POOL_SIZE;
2281  double *tempDifferent = new double[maxPool];
2282#define HASH 2
2283#if HASH > 1
2284  // for hashing
2285  typedef struct {
2286    int index, next;
2287  } CoinHashLink;
2288  CoinHashLink *hashThis = new CoinHashLink[4 * maxPool];
2289  for (int i = 0; i < 4 * maxPool; i++) {
2290    hashThis[i].index = -1;
2291    hashThis[i].next = -1;
2292  }
2293#endif
2294  int numberElements = 0;
2295  int hashDifferent = 0;
2296  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2297    CoinBigIndex k;
2298    columnStart_[iColumn] = numberElements;
2299    for (k = columnStart[iColumn]; k < columnStart[iColumn] + columnLength[iColumn];
2300         k++) {
2301      int iRow = row[k];
2302      double value = elementByColumn[k];
2303      poolInfo stuff;
2304      stuff.row_ = iRow;
2305#if HASH == 1 || HASH == 3
2306      int j;
2307      for (j = 0; j < numberDifferent_; j++) {
2308        if (value == tempDifferent[j])
2309          break;
2310      }
2311#endif
2312#if HASH == 2 || HASH == 3
2313      int ipos = hashit(value);
2314      int j1;
2315      while (true) {
2316        j1 = hashThis[ipos].index;
2317        if (j1 == -1) {
2318          hashThis[ipos].index = numberDifferent_;
2319          j1 = numberDifferent_;
2320          break;
2321        } else if (value == tempDifferent[j1]) {
2322          break;
2323        } else {
2324          int k = hashThis[ipos].next;
2325          if (k == -1) {
2326            j1 = numberDifferent_;
2327            while (true) {
2328              ++hashDifferent;
2329              if (hashThis[hashDifferent].index == -1) {
2330                break;
2331              }
2332            }
2333            hashThis[ipos].next = hashDifferent;
2334            hashThis[hashDifferent].index = numberDifferent_;
2335            break;
2336          } else {
2337            ipos = k;
2338            /* nothing worked - try it again */
2339          }
2340        }
2341      }
2342#if HASH == 2
2343      int j = j1;
2344#endif
2345#endif
2346#if HASH == 3
2347      assert(j == j1);
2348#endif
2349      if (j == numberDifferent_) {
2350        if (j == maxPool)
2351          break; // too many
2352        tempDifferent[j] = value;
2353        numberDifferent_++;
2354      }
2355      stuff.pool_ = j;
2356      stuff_[numberElements++] = stuff;
2357    }
2358  }
2359  columnStart_[numberColumns_] = numberElements;
2360#if HASH > 1
2361  delete[] hashThis;
2362#endif
2363  elements_ = new double[numberDifferent_];
2364  memcpy(elements_, tempDifferent, numberDifferent_ * sizeof(double));
2365  delete[] tempDifferent;
2366  if (numberDifferent_ == maxPool) {
2367    delete[] stuff_;
2368    delete[] elements_;
2369    delete[] columnStart_;
2370    elements_ = NULL;
2371    columnStart_ = NULL;
2372    stuff_ = NULL;
2373    numberRows_ = -1;
2374    numberColumns_ = -1;
2375    numberDifferent_ = -1;
2376  }
2377}
2378
2379//-------------------------------------------------------------------
2380// Destructor
2381//-------------------------------------------------------------------
2382ClpPoolMatrix::~ClpPoolMatrix()
2383{
2384  delete matrix_;
2385  delete[] elements_;
2386  delete[] columnStart_;
2387  delete[] lengths_;
2388  delete[] stuff_;
2389}
2390
2391//----------------------------------------------------------------
2392// Assignment operator
2393//-------------------------------------------------------------------
2394ClpPoolMatrix &
2395ClpPoolMatrix::operator=(const ClpPoolMatrix &rhs)
2396{
2397  if (this != &rhs) {
2398    ClpMatrixBase::operator=(rhs);
2399    delete matrix_;
2400    delete[] elements_;
2401    delete[] columnStart_;
2402    delete[] lengths_;
2403    delete[] stuff_;
2404    matrix_ = NULL;
2405    lengths_ = NULL;
2406    numberRows_ = rhs.numberRows_;
2407    numberColumns_ = rhs.numberColumns_;
2408    if (numberColumns_) {
2409      columnStart_ = CoinCopyOfArray(rhs.columnStart_, numberColumns_ + 1);
2410      CoinBigIndex numberElements = columnStart_[numberColumns_];
2411      stuff_ = CoinCopyOfArray(rhs.stuff_, numberElements);
2412      elements_ = CoinCopyOfArray(rhs.elements_, numberDifferent_);
2413    }
2414  }
2415  return *this;
2416}
2417// Constructor from arrays - handing over ownership
2418ClpPoolMatrix::ClpPoolMatrix(int numberColumns, CoinBigIndex *columnStart,
2419  poolInfo *stuff, double *elements)
2420  : ClpMatrixBase()
2421{
2422  setType(13);
2423  matrix_ = NULL;
2424  lengths_ = NULL;
2425  numberColumns_ = numberColumns;
2426  numberDifferent_ = 0;
2427  numberRows_ = 0;
2428  columnStart_ = columnStart;
2429  stuff_ = stuff;
2430  elements_ = elements;
2431  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2432    CoinBigIndex k;
2433    for (k = columnStart_[iColumn]; k < columnStart_[iColumn + 1];
2434         k++) {
2435      int iRow = stuff_[k].row_;
2436      int iPool = stuff_[k].pool_;
2437      numberDifferent_ = CoinMax(numberDifferent_, iPool);
2438      numberRows_ = CoinMax(numberRows_, iRow);
2439    }
2440  }
2441  // adjust
2442  numberDifferent_++;
2443  numberRows_++;
2444}
2445//-------------------------------------------------------------------
2446// Clone
2447//-------------------------------------------------------------------
2448ClpMatrixBase *ClpPoolMatrix::clone() const
2449{
2450  return new ClpPoolMatrix(*this);
2451}
2452/* Subset clone (without gaps).  Duplicates are allowed
2453   and order is as given */
2454ClpMatrixBase *
2455ClpPoolMatrix::subsetClone(int numberRows, const int *whichRows,
2456  int numberColumns,
2457  const int *whichColumns) const
2458{
2459  return new ClpPoolMatrix(*this, numberRows, whichRows,
2460    numberColumns, whichColumns);
2461}
2462/* Subset constructor (without gaps).  Duplicates are allowed
2463   and order is as given */
2464ClpPoolMatrix::ClpPoolMatrix(
2465  const ClpPoolMatrix &rhs,
2466  int numberRows, const int *whichRow,
2467  int numberColumns, const int *whichColumn)
2468  : ClpMatrixBase(rhs)
2469{
2470  setType(13);
2471  matrix_ = NULL;
2472  lengths_ = NULL;
2473  elements_ = NULL;
2474  columnStart_ = NULL;
2475  stuff_ = NULL;
2476  numberRows_ = numberRows;
2477  numberColumns_ = numberColumns;
2478  numberDifferent_ = 0;
2479  if (!numberRows_ || !numberColumns_)
2480    return;
2481  numberDifferent_ = rhs.numberDifferent_;
2482  elements_ = CoinCopyOfArray(rhs.elements_, numberDifferent_);
2483  int numberRowsOther = rhs.numberRows_;
2484  int numberColumnsOther = rhs.numberColumns_;
2485  int *newRow = new int[numberRowsOther + numberRows_];
2486  int *duplicateRow = newRow + numberRowsOther;
2487  for (int iRow = 0; iRow < numberRowsOther; iRow++)
2488    newRow[iRow] = -1;
2489  // and array for duplicating rows
2490  int numberBad = 0;
2491  for (int iRow = 0; iRow < numberRows_; iRow++) {
2492    duplicateRow[iRow] = -1;
2493    int kRow = whichRow[iRow];
2494    if (kRow >= 0 && kRow < numberRowsOther) {
2495      if (newRow[kRow] < 0) {
2496        // first time
2497        newRow[kRow] = iRow;
2498      } else {
2499        // duplicate
2500        int lastRow = newRow[kRow];
2501        newRow[kRow] = iRow;
2502        duplicateRow[iRow] = lastRow;
2503      }
2504    } else {
2505      // bad row
2506      numberBad++;
2507    }
2508  }
2509
2510  if (numberBad)
2511    throw CoinError("bad row entries",
2512      "subset constructor", "ClpPoolMatrix");
2513  // now get size and check columns
2514  CoinBigIndex size = 0;
2515  numberBad = 0;
2516  const CoinBigIndex *starts = rhs.columnStart_;
2517  const poolInfo *stuff = rhs.stuff_;
2518  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2519    int kColumn = whichColumn[iColumn];
2520    if (kColumn >= 0 && kColumn < numberColumnsOther) {
2521      CoinBigIndex i;
2522      for (i = starts[kColumn]; i < starts[kColumn + 1]; i++) {
2523        int kRow = stuff[i].row_;
2524        kRow = newRow[kRow];
2525        while (kRow >= 0) {
2526          size++;
2527          kRow = duplicateRow[kRow];
2528        }
2529      }
2530    } else {
2531      // bad column
2532      numberBad++;
2533    }
2534  }
2535  if (numberBad)
2536    throw CoinError("bad major entries",
2537      "subset constructor", "ClpPoolMatrix");
2538  // now create arrays
2539  stuff_ = new poolInfo[size];
2540  columnStart_ = new CoinBigIndex[numberColumns_ + 1];
2541  size = 0;
2542  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2543    int kColumn = whichColumn[iColumn];
2544    columnStart_[iColumn] = size;
2545    CoinBigIndex i;
2546    for (i = starts[kColumn]; i < starts[kColumn + 1]; i++) {
2547      int kRow = stuff[i].row_;
2548      int iPool = stuff[i].pool_;
2549      kRow = newRow[kRow];
2550      while (kRow >= 0) {
2551        stuff_[size].row_ = kRow;
2552        stuff_[size++].pool_ = iPool;
2553        kRow = duplicateRow[kRow];
2554      }
2555    }
2556  }
2557  columnStart_[numberColumns_] = size;
2558  delete[] newRow;
2559}
2560// Create matrix_
2561ClpPackedMatrix *
2562ClpPoolMatrix::createMatrix() const
2563{
2564  if (!matrix_) {
2565    CoinBigIndex numberElements = columnStart_[numberColumns_];
2566    double *elements = new double[numberElements];
2567    int *rows = new int[numberElements];
2568    bool lengthsExist = lengths_ != NULL;
2569    numberElements = 0;
2570    if (!lengthsExist)
2571      lengths_ = new int[numberColumns_];
2572    for (int i = 0; i < numberColumns_; i++) {
2573      lengths_[i] = static_cast< int >(columnStart_[i + 1] - columnStart_[i]);
2574      for (CoinBigIndex j = columnStart_[i]; j < columnStart_[i + 1]; j++) {
2575        elements[numberElements] = elements_[stuff_[j].pool_];
2576        rows[numberElements++] = stuff_[j].row_;
2577      }
2578    }
2579    CoinPackedMatrix *matrix = new CoinPackedMatrix(true, numberRows_, numberColumns_,
2580      numberElements,
2581      elements, rows,
2582      columnStart_, lengths_);
2583    delete[] elements;
2584    delete[] rows;
2585    matrix_ = new ClpPackedMatrix(matrix);
2586  }
2587  return matrix_;
2588}
2589
2590/* Returns a new matrix in reverse order without gaps */
2591ClpMatrixBase *
2592ClpPoolMatrix::reverseOrderedCopy() const
2593{
2594  // define to say no row copy
2595#define BIG_MATRIX
2596#ifndef BIG_MATRIX
2597  printf("XcreateMatrix at file %s line %d\n", __FILE__, __LINE__);
2598  createMatrix();
2599  return matrix_->reverseOrderedCopy();
2600#else
2601  return NULL;
2602#endif
2603}
2604#undef ABOCA_LITE
2605//unscaled versions
2606void ClpPoolMatrix::times(double scalar,
2607  const double *COIN_RESTRICT x, double *COIN_RESTRICT y) const
2608{
2609#if 0
2610  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2611  createMatrix()->times(scalar,x,y);
2612#else
2613  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2614    CoinBigIndex j;
2615    double value = x[iColumn];
2616    if (value) {
2617      CoinBigIndex start = columnStart_[iColumn];
2618      CoinBigIndex end = columnStart_[iColumn + 1];
2619      value *= scalar;
2620      for (j = start; j < end; j++) {
2621        int iRow = stuff_[j].row_;
2622        y[iRow] += value * elements_[stuff_[j].pool_];
2623      }
2624    }
2625  }
2626#endif
2627}
2628void ClpPoolMatrix::transposeTimes(double scalar,
2629  const double *COIN_RESTRICT x, double *COIN_RESTRICT y) const
2630{
2631#if 0
2632  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2633  createMatrix()->transposeTimes(scalar,x,y);
2634#else
2635  int iColumn;
2636  CoinBigIndex start = columnStart_[0];
2637  if (scalar == -1.0) {
2638    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2639      CoinBigIndex j;
2640      CoinBigIndex next = columnStart_[iColumn + 1];
2641      double value = 0.0;
2642      // scaled
2643      for (j = start; j < next; j++) {
2644        int jRow = stuff_[j].row_;
2645        value += x[jRow] * elements_[stuff_[j].pool_];
2646      }
2647      start = next;
2648      y[iColumn] -= value;
2649    }
2650  } else {
2651    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2652      CoinBigIndex j;
2653      CoinBigIndex next = columnStart_[iColumn + 1];
2654      double value = 0.0;
2655      // scaled
2656      for (j = start; j < next; j++) {
2657        int jRow = stuff_[j].row_;
2658        value += x[jRow] * elements_[stuff_[j].pool_];
2659      }
2660      start = next;
2661      y[iColumn] += value * scalar;
2662    }
2663  }
2664#endif
2665}
2666void ClpPoolMatrix::times(double scalar,
2667  const double *COIN_RESTRICT x, double *COIN_RESTRICT y,
2668  const double *rowScale,
2669  const double *columnScale) const
2670{
2671  printf("createMatrix at file %s line %d\n", __FILE__, __LINE__);
2672  createMatrix()->times(scalar, x, y);
2673}
2674void ClpPoolMatrix::transposeTimes(double scalar,
2675  const double *COIN_RESTRICT x, double *COIN_RESTRICT y,
2676  const double *rowScale,
2677  const double *columnScale,
2678  double *spare) const
2679{
2680#if 0
2681  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2682  createMatrix()->transposeTimes(scalar,x,y,rowScale,columnScale,spare);
2683#else
2684  if (rowScale) {
2685    int iColumn;
2686    if (!spare) {
2687      CoinBigIndex start = columnStart_[0];
2688      if (scalar == -1.0) {
2689        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2690          CoinBigIndex j;
2691          CoinBigIndex next = columnStart_[iColumn + 1];
2692          double value = 0.0;
2693          // scaled
2694          for (j = start; j < next; j++) {
2695            int jRow = stuff_[j].row_;
2696            value += x[jRow] * elements_[stuff_[j].pool_] * rowScale[jRow];
2697          }
2698          start = next;
2699          y[iColumn] -= value * columnScale[iColumn];
2700        }
2701      } else {
2702        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2703          CoinBigIndex j;
2704          CoinBigIndex next = columnStart_[iColumn + 1];
2705          double value = 0.0;
2706          // scaled
2707          for (j = start; j < next; j++) {
2708            int jRow = stuff_[j].row_;
2709            value += x[jRow] * elements_[stuff_[j].pool_] * rowScale[jRow];
2710          }
2711          start = next;
2712          y[iColumn] += value * scalar * columnScale[iColumn];
2713        }
2714      }
2715    } else {
2716      // can use spare region
2717      int iRow;
2718      int numberRows = matrix_->getNumRows();
2719      for (iRow = 0; iRow < numberRows; iRow++) {
2720        double value = x[iRow];
2721        if (value)
2722          spare[iRow] = value * rowScale[iRow];
2723        else
2724          spare[iRow] = 0.0;
2725      }
2726      CoinBigIndex start = columnStart_[0];
2727      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2728        CoinBigIndex j;
2729        CoinBigIndex next = columnStart_[iColumn + 1];
2730        double value = 0.0;
2731        // scaled
2732        for (j = start; j < next; j++) {
2733          int jRow = stuff_[j].row_;
2734          value += spare[jRow] * elements_[stuff_[j].pool_];
2735        }
2736        start = next;
2737        y[iColumn] += value * scalar * columnScale[iColumn];
2738      }
2739    }
2740  } else {
2741    transposeTimes(scalar, x, y);
2742  }
2743#endif
2744}
2745/* Return <code>x * A + y</code> in <code>z</code>.
2746   Squashes small elements and knows about ClpSimplex */
2747void ClpPoolMatrix::transposeTimes(const ClpSimplex *model, double scalar,
2748  const CoinIndexedVector *rowArray,
2749  CoinIndexedVector *y,
2750  CoinIndexedVector *columnArray) const
2751{
2752#if 0
2753  bool packed = rowArray->packedMode();
2754  assert (packed);
2755  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2756  createMatrix()->transposeTimes(model,scalar,rowArray,y,columnArray);
2757#else
2758  columnArray->clear();
2759  // do by column
2760  double *COIN_RESTRICT pi = rowArray->denseVector();
2761  int numberNonZero = 0;
2762  int *COIN_RESTRICT index = columnArray->getIndices();
2763  double *COIN_RESTRICT array = columnArray->denseVector();
2764  int numberInRowArray = rowArray->getNumElements();
2765  // maybe I need one in OsiSimplex
2766  double zeroTolerance = model->zeroTolerance();
2767  bool packed = rowArray->packedMode();
2768  assert(packed);
2769  // do by column
2770  const double *COIN_RESTRICT rowScale = model->rowScale();
2771  assert(!y->getNumElements());
2772  // need to expand pi into y
2773  assert(y->capacity() >= model->numberRows());
2774  double *COIN_RESTRICT piOld = pi;
2775  pi = y->denseVector();
2776  const int *COIN_RESTRICT whichRow = rowArray->getIndices();
2777  int i;
2778  // later combine phases
2779  if (!rowScale) {
2780    // modify pi so can collapse to one loop
2781    for (i = 0; i < numberInRowArray; i++) {
2782      int iRow = whichRow[i];
2783      pi[iRow] = scalar * piOld[i];
2784    }
2785    double value = 0.0;
2786    CoinBigIndex j;
2787    CoinBigIndex end = columnStart_[1];
2788    for (j = columnStart_[0]; j < end; j++) {
2789      int iRow = stuff_[j].row_;
2790      value += pi[iRow] * elements_[stuff_[j].pool_];
2791    }
2792    int iColumn;
2793    for (iColumn = 0; iColumn < numberColumns_ - 1; iColumn++) {
2794      CoinBigIndex start = end;
2795      end = columnStart_[iColumn + 2];
2796      if (fabs(value) > zeroTolerance) {
2797        array[numberNonZero] = value;
2798        index[numberNonZero++] = iColumn;
2799      }
2800      value = 0.0;
2801      for (j = start; j < end; j++) {
2802        int iRow = stuff_[j].row_;
2803        value += pi[iRow] * elements_[stuff_[j].pool_];
2804      }
2805    }
2806    if (fabs(value) > zeroTolerance) {
2807      array[numberNonZero] = value;
2808      index[numberNonZero++] = iColumn;
2809    }
2810  } else {
2811    // scaled
2812    // modify pi so can collapse to one loop
2813    for (i = 0; i < numberInRowArray; i++) {
2814      int iRow = whichRow[i];
2815      pi[iRow] = scalar * piOld[i] * rowScale[iRow];
2816    }
2817    const double *columnScale = model->columnScale();
2818    double value = 0.0;
2819    double scale = columnScale[0];
2820    CoinBigIndex j;
2821    CoinBigIndex end = columnStart_[1];
2822    for (j = columnStart_[0]; j < end; j++) {
2823      int iRow = stuff_[j].row_;
2824      value += pi[iRow] * elements_[stuff_[j].pool_];
2825    }
2826    int iColumn;
2827    for (iColumn = 0; iColumn < numberColumns_ - 1; iColumn++) {
2828      value *= scale;
2829      CoinBigIndex start = end;
2830      scale = columnScale[iColumn + 1];
2831      end = columnStart_[iColumn + 2];
2832      if (fabs(value) > zeroTolerance) {
2833        array[numberNonZero] = value;
2834        index[numberNonZero++] = iColumn;
2835      }
2836      value = 0.0;
2837      for (j = start; j < end; j++) {
2838        int iRow = stuff_[j].row_;
2839        value += pi[iRow] * elements_[stuff_[j].pool_];
2840      }
2841    }
2842    value *= scale;
2843    if (fabs(value) > zeroTolerance) {
2844      array[numberNonZero] = value;
2845      index[numberNonZero++] = iColumn;
2846    }
2847  }
2848  // zero out
2849  int numberRows = model->numberRows();
2850  if (numberInRowArray * 4 < numberRows) {
2851    for (i = 0; i < numberInRowArray; i++) {
2852      int iRow = whichRow[i];
2853      pi[iRow] = 0.0;
2854    }
2855  } else {
2856    CoinZeroN(pi, numberRows);
2857  }
2858  columnArray->setNumElements(numberNonZero);
2859  y->setNumElements(0);
2860  columnArray->setPackedMode(true);
2861#endif
2862}
2863/* Return <code>x * A + y</code> in <code>z</code>.
2864   Squashes small elements and knows about ClpSimplex */
2865void ClpPoolMatrix::transposeTimesByRow(const ClpSimplex *model, double scalar,
2866  const CoinIndexedVector *rowArray,
2867  CoinIndexedVector *y,
2868  CoinIndexedVector *columnArray) const
2869{
2870  printf("createMatrix at file %s line %d\n", __FILE__, __LINE__);
2871  createMatrix()->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
2872}
2873/* Return <code>x *A in <code>z</code> but
2874   just for indices in y. */
2875void ClpPoolMatrix::subsetTransposeTimes(const ClpSimplex *model,
2876  const CoinIndexedVector *rowArray,
2877  const CoinIndexedVector *y,
2878  CoinIndexedVector *columnArray) const
2879{
2880#if 0
2881  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2882  createMatrix()->subsetTransposeTimes(model,rowArray,y,columnArray);
2883#else
2884  columnArray->clear();
2885  double *COIN_RESTRICT pi = rowArray->denseVector();
2886  double *COIN_RESTRICT array = columnArray->denseVector();
2887  int jColumn;
2888  // get matrix data pointers
2889  const double *COIN_RESTRICT rowScale = model->rowScale();
2890  int numberToDo = y->getNumElements();
2891  const int *COIN_RESTRICT which = y->getIndices();
2892  assert(!rowArray->packedMode());
2893  columnArray->setPacked();
2894  if (!rowScale) {
2895    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
2896      int iColumn = which[jColumn];
2897      double value = 0.0;
2898      CoinBigIndex start = columnStart_[iColumn];
2899      CoinBigIndex end = columnStart_[iColumn + 1];
2900      array[jColumn] = value;
2901      CoinBigIndex j;
2902      for (j = start; j < end; j++) {
2903        int iRow = stuff_[j].row_;
2904        value += pi[iRow] * elements_[stuff_[j].pool_];
2905      }
2906      array[jColumn] = value;
2907    }
2908  } else {
2909    // scaled
2910    const double *columnScale = model->columnScale();
2911    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
2912      int iColumn = which[jColumn];
2913      double value = 0.0;
2914      CoinBigIndex start = columnStart_[iColumn];
2915      CoinBigIndex end = columnStart_[iColumn + 1];
2916      CoinBigIndex j;
2917      for (j = start; j < end; j++) {
2918        int iRow = stuff_[j].row_;
2919        value += pi[iRow] * elements_[stuff_[j].pool_] * rowScale[iRow];
2920      }
2921      array[jColumn] = value * columnScale[iColumn];
2922    }
2923  }
2924#endif
2925}
2926/// returns number of elements in column part of basis,
2927int ClpPoolMatrix::countBasis(const int *whichColumn,
2928  int &numberColumnBasic)
2929{
2930  int i;
2931  CoinBigIndex numberElements = 0;
2932  for (i = 0; i < numberColumnBasic; i++) {
2933    int iColumn = whichColumn[i];
2934    numberElements += columnStart_[iColumn + 1] - columnStart_[iColumn];
2935  }
2936#if COIN_BIG_INDEX
2937  if (numberElements > COIN_INT_MAX) {
2938    printf("Factorization too large\n");
2939    abort();
2940  }
2941#endif
2942  return static_cast< int >(numberElements);
2943}
2944void ClpPoolMatrix::fillBasis(ClpSimplex *,
2945  const int *whichColumn,
2946  int &numberColumnBasic,
2947  int *indexRowU, int *start,
2948  int *rowCount, int *columnCount,
2949  CoinFactorizationDouble *elementU)
2950{
2951  CoinBigIndex numberElements = start[0];
2952  for (int i = 0; i < numberColumnBasic; i++) {
2953    int iColumn = whichColumn[i];
2954    CoinBigIndex j = columnStart_[iColumn];
2955    for (; j < columnStart_[iColumn + 1]; j++) {
2956      int iRow = stuff_[j].row_;
2957      indexRowU[numberElements] = iRow;
2958      rowCount[iRow]++;
2959      elementU[numberElements++] = elements_[stuff_[j].pool_];
2960      ;
2961    }
2962    start[i + 1] = static_cast< int >(numberElements);
2963    columnCount[i] = static_cast< int >(numberElements - start[i]);
2964  }
2965#if COIN_BIG_INDEX
2966  if (numberElements > COIN_INT_MAX)
2967    abort();
2968#endif
2969}
2970/* Unpacks a column into an CoinIndexedvector
2971 */
2972void ClpPoolMatrix::unpack(const ClpSimplex *,
2973  CoinIndexedVector *rowArray,
2974  int iColumn) const
2975{
2976  CoinBigIndex j = columnStart_[iColumn];
2977  for (; j < columnStart_[iColumn + 1]; j++) {
2978    int iRow = stuff_[j].row_;
2979    rowArray->add(iRow, elements_[stuff_[j].pool_]);
2980  }
2981}
2982/* Unpacks a column into an CoinIndexedvector
2983** in packed foramt
2984Note that model is NOT const.  Bounds and objective could
2985be modified if doing column generation (just for this variable) */
2986void ClpPoolMatrix::unpackPacked(ClpSimplex *,
2987  CoinIndexedVector *rowArray,
2988  int iColumn) const
2989{
2990  int *COIN_RESTRICT index = rowArray->getIndices();
2991  double *COIN_RESTRICT array = rowArray->denseVector();
2992  int number = 0;
2993  CoinBigIndex j = columnStart_[iColumn];
2994  for (; j < columnStart_[iColumn + 1]; j++) {
2995    int iRow = stuff_[j].row_;
2996    array[number] = elements_[stuff_[j].pool_];
2997    index[number++] = iRow;
2998  }
2999  rowArray->setNumElements(number);
3000  rowArray->setPackedMode(true);
3001}
3002/* Adds multiple of a column into an CoinIndexedvector
3003   You can use quickAdd to add to vector */
3004void ClpPoolMatrix::add(const ClpSimplex *, CoinIndexedVector *rowArray,
3005  int iColumn, double multiplier) const
3006{
3007  CoinBigIndex j = columnStart_[iColumn];
3008  for (; j < columnStart_[iColumn + 1]; j++) {
3009    int iRow = stuff_[j].row_;
3010    rowArray->quickAdd(iRow, multiplier * elements_[stuff_[j].pool_]);
3011  }
3012}
3013/* Adds multiple of a column into an array */
3014void ClpPoolMatrix::add(const ClpSimplex *, double *array,
3015  int iColumn, double multiplier) const
3016{
3017  CoinBigIndex j = columnStart_[iColumn];
3018  for (; j < columnStart_[iColumn + 1]; j++) {
3019    int iRow = stuff_[j].row_;
3020    array[iRow] += multiplier * elements_[stuff_[j].pool_];
3021  }
3022}
3023
3024// Return a complete CoinPackedMatrix
3025CoinPackedMatrix *
3026ClpPoolMatrix::getPackedMatrix() const
3027{
3028  createMatrix();
3029  return matrix_->matrix();
3030}
3031/* A vector containing the elements in the packed matrix. Note that there
3032   might be gaps in this list, entries that do not belong to any
3033   major-dimension vector. To get the actual elements one should look at
3034   this vector together with vectorStarts and vectorLengths. */
3035const double *
3036ClpPoolMatrix::getElements() const
3037{
3038  printf("createMatrix at file %s line %d\n", __FILE__, __LINE__);
3039  return createMatrix()->getElements();
3040}
3041
3042/* A vector containing the minor indices of the elements in the packed
3043   matrix. Note that there might be gaps in this list, entries that do not
3044   belong to any major-dimension vector. To get the actual elements one
3045   should look at this vector together with vectorStarts and
3046   vectorLengths. */
3047const int *
3048ClpPoolMatrix::getIndices() const
3049{
3050  printf("createMatrix at file %s line %d\n", __FILE__, __LINE__);
3051  return createMatrix()->getIndices();
3052}
3053const CoinBigIndex *
3054ClpPoolMatrix::getVectorStarts() const
3055{
3056  printf("createMatrix at file %s line %d\n", __FILE__, __LINE__);
3057  return createMatrix()->getVectorStarts();
3058}
3059/* The lengths of the major-dimension vectors. */
3060const int *
3061ClpPoolMatrix::getVectorLengths() const
3062{
3063#if 0
3064  printf("XcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
3065  return createMatrix()->getVectorLengths();
3066#else
3067  if (!lengths_) {
3068    lengths_ = new int[numberColumns_];
3069    for (int i = 0; i < numberColumns_; i++)
3070      lengths_[i] = static_cast< int >(columnStart_[i + 1] - columnStart_[i]);
3071  }
3072  return lengths_;
3073#endif
3074}
3075/* The length of a major-dimension vector. */
3076int ClpPoolMatrix::getVectorLength(int index) const
3077{
3078  return static_cast< int >(columnStart_[index + 1] - columnStart_[index]);
3079}
3080/* Delete the columns whose indices are listed in <code>indDel</code>. */
3081void ClpPoolMatrix::deleteCols(const int numDel, const int *indDel)
3082{
3083  int iColumn;
3084  CoinBigIndex newSize = columnStart_[numberColumns_];
3085  ;
3086  int numberBad = 0;
3087  // Use array to make sure we can have duplicates
3088  int *which = new int[numberColumns_];
3089  memset(which, 0, numberColumns_ * sizeof(int));
3090  int nDuplicate = 0;
3091  for (iColumn = 0; iColumn < numDel; iColumn++) {
3092    int jColumn = indDel[iColumn];
3093    if (jColumn < 0 || jColumn >= numberColumns_) {
3094      numberBad++;
3095    } else {
3096      newSize -= columnStart_[jColumn + 1] - columnStart_[jColumn];
3097      if (which[jColumn])
3098        nDuplicate++;
3099      else
3100        which[jColumn] = 1;
3101    }
3102  }
3103  if (numberBad)
3104    throw CoinError("Indices out of range", "deleteCols", "ClpPoolMatrix");
3105  int newNumber = numberColumns_ - numDel + nDuplicate;
3106  // Get rid of temporary arrays
3107  delete matrix_;
3108  matrix_ = NULL;
3109  CoinBigIndex *columnStart = new CoinBigIndex[newNumber + 1];
3110  poolInfo *stuff = new poolInfo[newSize];
3111  newNumber = 0;
3112  newSize = 0;
3113  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3114    if (!which[iColumn]) {
3115      CoinBigIndex start, end;
3116      CoinBigIndex i;
3117      start = columnStart_[iColumn];
3118      end = columnStart_[iColumn + 1];
3119      columnStart[newNumber] = newSize;
3120      for (i = start; i < end; i++)
3121        stuff[newSize++] = stuff_[i];
3122    }
3123  }
3124  columnStart[newNumber] = newSize;
3125  delete[] which;
3126  delete[] columnStart_;
3127  columnStart_ = columnStart;
3128  delete[] stuff_;
3129  stuff_ = stuff;
3130  numberColumns_ = newNumber;
3131}
3132/* Delete the rows whose indices are listed in <code>indDel</code>. */
3133void ClpPoolMatrix::deleteRows(const int numDel, const int *indDel)
3134{
3135  int iRow;
3136  int numberBad = 0;
3137  // Use array to make sure we can have duplicates
3138  int *which = new int[numberRows_];
3139  memset(which, 0, numberRows_ * sizeof(int));
3140  int nDuplicate = 0;
3141  for (iRow = 0; iRow < numDel; iRow++) {
3142    int jRow = indDel[iRow];
3143    if (jRow < 0 || jRow >= numberRows_) {
3144      numberBad++;
3145    } else {
3146      if (which[jRow])
3147        nDuplicate++;
3148      else
3149        which[jRow] = 1;
3150    }
3151  }
3152  if (numberBad)
3153    throw CoinError("Indices out of range", "deleteRows", "ClpPoolMatrix");
3154  CoinBigIndex iElement;
3155  CoinBigIndex numberElements = columnStart_[numberColumns_];
3156  CoinBigIndex newSize = 0;
3157  for (iElement = 0; iElement < numberElements; iElement++) {
3158    iRow = stuff_[iElement].row_;
3159    if (!which[iRow])
3160      newSize++;
3161  }
3162  int newNumber = numberRows_ - numDel + nDuplicate;
3163  // Get rid of temporary arrays
3164  delete matrix_;
3165  matrix_ = NULL;
3166  poolInfo *stuff = new poolInfo[newSize];
3167  newSize = 0;
3168  // redo which
3169  int numberRows = 0;
3170  for (iRow = 0; iRow < numberRows_; iRow++) {
3171    if (which[iRow]) {
3172      which[iRow] = -1; // not wanted
3173    } else {
3174      which[iRow] = numberRows;
3175      numberRows++;
3176    }
3177  }
3178  int iColumn;
3179  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3180    CoinBigIndex start, end;
3181    CoinBigIndex i;
3182    start = columnStart_[iColumn];
3183    end = columnStart_[iColumn + 1];
3184    columnStart_[newNumber] = newSize;
3185    for (i = start; i < end; i++) {
3186      poolInfo newStuff = stuff_[i];
3187      iRow = which[newStuff.row_];
3188      if (iRow >= 0) {
3189        newStuff.row_ = iRow;
3190        stuff[newSize++] = newStuff;
3191      }
3192    }
3193  }
3194  columnStart_[numberColumns_] = newSize;
3195  delete[] which;
3196  delete[] stuff_;
3197  stuff_ = stuff;
3198  numberRows_ = newNumber;
3199}
3200bool ClpPoolMatrix::isColOrdered() const
3201{
3202  return true;
3203}
3204/* Number of entries in the packed matrix. */
3205CoinBigIndex
3206ClpPoolMatrix::getNumElements() const
3207{
3208  if (numberColumns_)
3209    return columnStart_[numberColumns_];
3210  else
3211    return 0;
3212}
3213/* Returns largest and smallest elements of both signs.
3214   Largest refers to largest absolute value.
3215*/
3216void ClpPoolMatrix::rangeOfElements(double &smallestNegative, double &largestNegative,
3217  double &smallestPositive, double &largestPositive)
3218{
3219  smallestNegative = -COIN_DBL_MAX;
3220  largestNegative = 0.0;
3221  smallestPositive = COIN_DBL_MAX;
3222  largestPositive = 0.0;
3223  for (int i = 0; i < numberDifferent_; i++) {
3224    double value = elements_[i];
3225    if (value > 0.0) {
3226      smallestPositive = CoinMin(smallestPositive, value);
3227      largestPositive = CoinMax(largestPositive, value);
3228    } else if (value < 0.0) {
3229      smallestNegative = CoinMax(smallestNegative, value);
3230      largestNegative = CoinMin(largestNegative, value);
3231    }
3232  }
3233}
3234// Says whether it can do partial pricing
3235bool ClpPoolMatrix::canDoPartialPricing() const
3236{
3237  return true;
3238}
3239// Partial pricing
3240void ClpPoolMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
3241  int &bestSequence, int &numberWanted)
3242{
3243  printf("createMatrix at file %s line %d\n", __FILE__, __LINE__);
3244  createMatrix()->partialPricing(model, startFraction,
3245    endFraction, bestSequence,
3246    numberWanted);
3247}
3248// Allow any parts of a created CoinMatrix to be deleted
3249void ClpPoolMatrix::releasePackedMatrix() const
3250{
3251  delete matrix_;
3252  matrix_ = NULL;
3253}
3254// These have to match ClpPrimalColumnSteepest version
3255#define reference(i) (((reference[i >> 5] >> (i & 31)) & 1) != 0)
3256/* Updates two arrays for steepest and does devex weights
3257   Returns nonzero if updates reduced cost and infeas -
3258   new infeas in dj1  - This does not at present*/
3259int ClpPoolMatrix::transposeTimes2(const ClpSimplex *model,
3260  const CoinIndexedVector *pi1, CoinIndexedVector *dj1,
3261  const CoinIndexedVector *pi2,
3262  CoinIndexedVector *spare,
3263  double *infeas, double *reducedCost,
3264  double referenceIn, double devex,
3265  // Array for exact devex to say what is in reference framework
3266  unsigned int *COIN_RESTRICT reference,
3267  double *COIN_RESTRICT weights, double scaleFactor)
3268{
3269#if 0
3270  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
3271  return createMatrix()->transposeTimes2(model,pi1,dj1,pi2,spare,
3272                                         infeas,reducedCost,referenceIn,devex,
3273                                         reference,weights,scaleFactor);
3274#else
3275  int returnCode = 0;
3276  // put row of tableau in dj1
3277  double *COIN_RESTRICT pi = pi1->denseVector();
3278  int numberNonZero = 0;
3279  int *COIN_RESTRICT index = dj1->getIndices();
3280  double *COIN_RESTRICT array = dj1->denseVector();
3281  int numberInRowArray = pi1->getNumElements();
3282  double zeroTolerance = model->zeroTolerance();
3283  double dualTolerance = model->currentDualTolerance();
3284  // we can't really trust infeasibilities if there is dual error
3285  // this coding has to mimic coding in checkDualSolution
3286  double error = CoinMin(1.0e-2, model->largestDualError());
3287  // allow tolerance at least slightly bigger than standard
3288  dualTolerance = dualTolerance + error;
3289  bool packed = pi1->packedMode();
3290  assert(packed);
3291  // do by column
3292  int iColumn;
3293  const double *COIN_RESTRICT rowScale = model->rowScale();
3294  assert(!spare->getNumElements());
3295  assert(numberColumns_ > 0);
3296  double *COIN_RESTRICT piWeight = pi2->denseVector();
3297  assert(!pi2->packedMode());
3298  bool killDjs = (scaleFactor == 0.0);
3299  if (!scaleFactor)
3300    scaleFactor = 1.0;
3301  // need to expand pi into y
3302  assert(spare->capacity() >= model->numberRows());
3303  double *COIN_RESTRICT piOld = pi;
3304  pi = spare->denseVector();
3305  const int *COIN_RESTRICT whichRow = pi1->getIndices();
3306  int i;
3307  if (!rowScale) {
3308    // modify pi so can collapse to one loop
3309    for (i = 0; i < numberInRowArray; i++) {
3310      int iRow = whichRow[i];
3311      pi[iRow] = piOld[i];
3312    }
3313    if (infeas)
3314      returnCode = 1;
3315    CoinBigIndex j;
3316    CoinBigIndex end = columnStart_[0];
3317    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3318      CoinBigIndex start = end;
3319      end = columnStart_[iColumn + 1];
3320      ClpSimplex::Status status = model->getStatus(iColumn);
3321      if (status == ClpSimplex::basic || status == ClpSimplex::isFixed)
3322        continue;
3323      double value = 0.0;
3324      for (j = start; j < end; j++) {
3325        int iRow = stuff_[j].row_;
3326        value -= pi[iRow] * elements_[stuff_[j].pool_];
3327      }
3328      if (fabs(value) > zeroTolerance) {
3329        // and do other array
3330        double modification = 0.0;
3331        for (j = start; j < end; j++) {
3332          int iRow = stuff_[j].row_;
3333          modification += piWeight[iRow] * elements_[stuff_[j].pool_];
3334        }
3335        double thisWeight = weights[iColumn];
3336        double pivot = value * scaleFactor;
3337        double pivotSquared = pivot * pivot;
3338        thisWeight += pivotSquared * devex + pivot * modification;
3339        if (thisWeight < DEVEX_TRY_NORM) {
3340          if (referenceIn < 0.0) {
3341            // steepest
3342            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
3343          } else {
3344            // exact
3345            thisWeight = referenceIn * pivotSquared;
3346            if (reference(iColumn))
3347              thisWeight += 1.0;
3348            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
3349          }
3350        }
3351        weights[iColumn] = thisWeight;
3352        if (!killDjs) {
3353          value = reducedCost[iColumn] - value;
3354          reducedCost[iColumn] = value;
3355          // simplify status
3356          ClpSimplex::Status status = model->getStatus(iColumn);
3357
3358          switch (status) {
3359
3360          case ClpSimplex::basic:
3361          case ClpSimplex::isFixed:
3362            break;
3363          case ClpSimplex::isFree:
3364          case ClpSimplex::superBasic:
3365            if (fabs(value) > FREE_ACCEPT * dualTolerance) {
3366              // we are going to bias towards free (but only if reasonable)
3367              value *= FREE_BIAS;
3368              value *= value;
3369              // store square in list
3370              if (infeas[iColumn]) {
3371                infeas[iColumn] = value; // already there
3372              } else {
3373                array[numberNonZero] = value;
3374                index[numberNonZero++] = iColumn;
3375              }
3376            } else {
3377              array[numberNonZero] = 0.0;
3378              index[numberNonZero++] = iColumn;
3379            }
3380            break;
3381          case ClpSimplex::atUpperBound:
3382            if (value > dualTolerance) {
3383              value *= value;
3384              // store square in list
3385              if (infeas[iColumn]) {
3386                infeas[iColumn] = value; // already there
3387              } else {
3388                array[numberNonZero] = value;
3389                index[numberNonZero++] = iColumn;
3390              }
3391            } else {
3392              array[numberNonZero] = 0.0;
3393              index[numberNonZero++] = iColumn;
3394            }
3395            break;
3396          case ClpSimplex::atLowerBound:
3397            if (value < -dualTolerance) {
3398              value *= value;
3399              // store square in list
3400              if (infeas[iColumn]) {
3401                infeas[iColumn] = value; // already there
3402              } else {
3403                array[numberNonZero] = value;
3404                index[numberNonZero++] = iColumn;
3405              }
3406            } else {
3407              array[numberNonZero] = 0.0;
3408              index[numberNonZero++] = iColumn;
3409            }
3410          }
3411        }
3412      }
3413    }
3414  } else {
3415    // scaled
3416    // modify pi so can collapse to one loop
3417    for (i = 0; i < numberInRowArray; i++) {
3418      int iRow = whichRow[i];
3419      pi[iRow] = piOld[i] * rowScale[iRow];
3420    }
3421    // can also scale piWeight as not used again
3422    int numberWeight = pi2->getNumElements();
3423    const int *indexWeight = pi2->getIndices();
3424    for (i = 0; i < numberWeight; i++) {
3425      int iRow = indexWeight[i];
3426      piWeight[iRow] *= rowScale[iRow];
3427    }
3428    if (infeas)
3429      returnCode = 1;
3430    const double *COIN_RESTRICT columnScale = model->columnScale();
3431    CoinBigIndex j;
3432    CoinBigIndex end = columnStart_[0];
3433    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3434      CoinBigIndex start = end;
3435      end = columnStart_[iColumn + 1];
3436      ClpSimplex::Status status = model->getStatus(iColumn);
3437      if (status == ClpSimplex::basic || status == ClpSimplex::isFixed)
3438        continue;
3439      double scale = columnScale[iColumn];
3440      double value = 0.0;
3441      for (j = start; j < end; j++) {
3442        int iRow = stuff_[j].row_;
3443        value -= pi[iRow] * elements_[stuff_[j].pool_];
3444      }
3445      value *= scale;
3446      if (fabs(value) > zeroTolerance) {
3447        double modification = 0.0;
3448        for (j = start; j < end; j++) {
3449          int iRow = stuff_[j].row_;
3450          modification += piWeight[iRow] * elements_[stuff_[j].pool_];
3451        }
3452        modification *= scale;
3453        double thisWeight = weights[iColumn];
3454        double pivot = value * scaleFactor;
3455        double pivotSquared = pivot * pivot;
3456        thisWeight += pivotSquared * devex + pivot * modification;
3457        if (thisWeight < DEVEX_TRY_NORM) {
3458          if (referenceIn < 0.0) {
3459            // steepest
3460            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
3461          } else {
3462            // exact
3463            thisWeight = referenceIn * pivotSquared;
3464            if (reference(iColumn))
3465              thisWeight += 1.0;
3466            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
3467          }
3468        }
3469        weights[iColumn] = thisWeight;
3470        if (!killDjs) {
3471          value = reducedCost[iColumn] - value;
3472          reducedCost[iColumn] = value;
3473          // simplify status
3474          ClpSimplex::Status status = model->getStatus(iColumn);
3475
3476          switch (status) {
3477
3478          case ClpSimplex::basic:
3479          case ClpSimplex::isFixed:
3480            break;
3481          case ClpSimplex::isFree:
3482          case ClpSimplex::superBasic:
3483            if (fabs(value) > FREE_ACCEPT * dualTolerance) {
3484              // we are going to bias towards free (but only if reasonable)
3485              value *= FREE_BIAS;
3486              value *= value;
3487              // store square in list
3488              if (infeas[iColumn]) {
3489                infeas[iColumn] = value; // already there
3490              } else {
3491                array[numberNonZero] = value;
3492                index[numberNonZero++] = iColumn;
3493              }
3494            } else {
3495              array[numberNonZero] = 0.0;
3496              index[numberNonZero++] = iColumn;
3497            }
3498            break;
3499          case ClpSimplex::atUpperBound:
3500            if (value > dualTolerance) {
3501              value *= value;
3502              // store square in list
3503              if (infeas[iColumn]) {
3504                infeas[iColumn] = value; // already there
3505              } else {
3506                array[numberNonZero] = value;
3507                index[numberNonZero++] = iColumn;
3508              }
3509            } else {
3510              array[numberNonZero] = 0.0;
3511              index[numberNonZero++] = iColumn;
3512            }
3513            break;
3514          case ClpSimplex::atLowerBound:
3515            if (value < -dualTolerance) {
3516              value *= value;
3517              // store square in list
3518              if (infeas[iColumn]) {
3519                infeas[iColumn] = value; // already there
3520              } else {
3521                array[numberNonZero] = value;
3522                index[numberNonZero++] = iColumn;
3523              }
3524            } else {
3525              array[numberNonZero] = 0.0;
3526              index[numberNonZero++] = iColumn;
3527            }
3528          }
3529        }
3530      }
3531    }
3532  }
3533  // zero out
3534  int numberRows = model->numberRows();
3535  if (numberInRowArray * 4 < numberRows) {
3536    for (i = 0; i < numberInRowArray; i++) {
3537      int iRow = whichRow[i];
3538      pi[iRow] = 0.0;
3539    }
3540  } else {
3541    CoinZeroN(pi, numberRows);
3542  }
3543  dj1->setNumElements(numberNonZero);
3544  spare->setNumElements(0);
3545  if (packed)
3546    dj1->setPackedMode(true);
3547  return returnCode;
3548#endif
3549}
3550// Updates second array for steepest and does devex weights
3551void ClpPoolMatrix::subsetTimes2(const ClpSimplex *model,
3552  CoinIndexedVector *dj1,
3553  const CoinIndexedVector *pi2, CoinIndexedVector *dj2,
3554  double referenceIn, double devex,
3555  // Array for exact devex to say what is in reference framework
3556  unsigned int *COIN_RESTRICT reference,
3557  double *COIN_RESTRICT weights, double scaleFactor)
3558{
3559  printf("createMatrix at file %s line %d\n", __FILE__, __LINE__);
3560  createMatrix()->subsetTimes2(model, dj1, pi2, dj2,
3561    referenceIn, devex,
3562    reference, weights, scaleFactor);
3563}
3564/* Set the dimensions of the matrix. In effect, append new empty
3565   columns/rows to the matrix. A negative number for either dimension
3566   means that that dimension doesn't change. Otherwise the new dimensions
3567   MUST be at least as large as the current ones otherwise an exception
3568   is thrown. */
3569void ClpPoolMatrix::setDimensions(int newnumrows, int newnumcols)
3570{
3571  // at present abort
3572  abort();
3573}
3574#endif
3575
3576/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
3577*/
Note: See TracBrowser for help on using the repository browser.