source: trunk/Clp/src/ClpNetworkMatrix.cpp

Last change on this file 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: 36.3 KB
Line 
1/* $Id: ClpNetworkMatrix.cpp 2385 2019-01-06 19:43:06Z forrest $ */
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 "CoinHelperFunctions.hpp"
11#include "CoinPackedVector.hpp"
12
13#include "ClpSimplex.hpp"
14#include "ClpFactorization.hpp"
15// at end to get min/max!
16#include "ClpNetworkMatrix.hpp"
17#include "ClpPlusMinusOneMatrix.hpp"
18#include "ClpMessage.hpp"
19#include <iostream>
20#include <cassert>
21
22//#############################################################################
23// Constructors / Destructor / Assignment
24//#############################################################################
25
26//-------------------------------------------------------------------
27// Default Constructor
28//-------------------------------------------------------------------
29ClpNetworkMatrix::ClpNetworkMatrix()
30  : ClpMatrixBase()
31{
32  setType(11);
33  matrix_ = NULL;
34  lengths_ = NULL;
35  indices_ = NULL;
36  numberRows_ = 0;
37  numberColumns_ = 0;
38  trueNetwork_ = false;
39}
40
41/* Constructor from two arrays */
42ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int *head,
43  const int *tail)
44  : ClpMatrixBase()
45{
46  setType(11);
47  matrix_ = NULL;
48  lengths_ = NULL;
49  indices_ = new int[2 * numberColumns];
50  ;
51  numberRows_ = -1;
52  numberColumns_ = numberColumns;
53  trueNetwork_ = true;
54  int iColumn;
55  CoinBigIndex j = 0;
56  for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
57    int iRow = head[iColumn];
58    numberRows_ = CoinMax(numberRows_, iRow);
59    indices_[j] = iRow;
60    iRow = tail[iColumn];
61    numberRows_ = CoinMax(numberRows_, iRow);
62    indices_[j + 1] = iRow;
63  }
64  numberRows_++;
65}
66//-------------------------------------------------------------------
67// Copy constructor
68//-------------------------------------------------------------------
69ClpNetworkMatrix::ClpNetworkMatrix(const ClpNetworkMatrix &rhs)
70  : ClpMatrixBase(rhs)
71{
72  matrix_ = NULL;
73  lengths_ = NULL;
74  indices_ = NULL;
75  numberRows_ = rhs.numberRows_;
76  numberColumns_ = rhs.numberColumns_;
77  trueNetwork_ = rhs.trueNetwork_;
78  if (numberColumns_) {
79    indices_ = new int[2 * numberColumns_];
80    CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
81  }
82  int numberRows = getNumRows();
83  if (rhs.rhsOffset_ && numberRows) {
84    rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
85  } else {
86    rhsOffset_ = NULL;
87  }
88}
89
90ClpNetworkMatrix::ClpNetworkMatrix(const CoinPackedMatrix &rhs)
91  : ClpMatrixBase()
92{
93  setType(11);
94  matrix_ = NULL;
95  lengths_ = NULL;
96  indices_ = NULL;
97  int iColumn;
98  assert(rhs.isColOrdered());
99  // get matrix data pointers
100  const int *row = rhs.getIndices();
101  const CoinBigIndex *columnStart = rhs.getVectorStarts();
102  const int *columnLength = rhs.getVectorLengths();
103  const double *elementByColumn = rhs.getElements();
104  numberColumns_ = rhs.getNumCols();
105  int goodNetwork = 1;
106  numberRows_ = -1;
107  indices_ = new int[2 * numberColumns_];
108  CoinBigIndex j = 0;
109  for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
110    CoinBigIndex k = columnStart[iColumn];
111    int iRow;
112    switch (columnLength[iColumn]) {
113    case 0:
114      goodNetwork = -1; // not classic network
115      indices_[j] = -1;
116      indices_[j + 1] = -1;
117      break;
118
119    case 1:
120      goodNetwork = -1; // not classic network
121      if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
122        indices_[j] = -1;
123        iRow = row[k];
124        numberRows_ = CoinMax(numberRows_, iRow);
125        indices_[j + 1] = iRow;
126      } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
127        indices_[j + 1] = -1;
128        iRow = row[k];
129        numberRows_ = CoinMax(numberRows_, iRow);
130        indices_[j] = iRow;
131      } else {
132        goodNetwork = 0; // not a network
133      }
134      break;
135
136    case 2:
137      if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
138        if (fabs(elementByColumn[k + 1] + 1.0) < 1.0e-10) {
139          iRow = row[k];
140          numberRows_ = CoinMax(numberRows_, iRow);
141          indices_[j + 1] = iRow;
142          iRow = row[k + 1];
143          numberRows_ = CoinMax(numberRows_, iRow);
144          indices_[j] = iRow;
145        } else {
146          goodNetwork = 0; // not a network
147        }
148      } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
149        if (fabs(elementByColumn[k + 1] - 1.0) < 1.0e-10) {
150          iRow = row[k];
151          numberRows_ = CoinMax(numberRows_, iRow);
152          indices_[j] = iRow;
153          iRow = row[k + 1];
154          numberRows_ = CoinMax(numberRows_, iRow);
155          indices_[j + 1] = iRow;
156        } else {
157          goodNetwork = 0; // not a network
158        }
159      } else {
160        goodNetwork = 0; // not a network
161      }
162      break;
163
164    default:
165      goodNetwork = 0; // not a network
166      break;
167    }
168    if (!goodNetwork)
169      break;
170  }
171  if (!goodNetwork) {
172    delete[] indices_;
173    // put in message
174    printf("Not a network - can test if indices_ null\n");
175    indices_ = NULL;
176    numberRows_ = 0;
177    numberColumns_ = 0;
178  } else {
179    numberRows_++; //  correct
180    trueNetwork_ = goodNetwork > 0;
181  }
182}
183
184//-------------------------------------------------------------------
185// Destructor
186//-------------------------------------------------------------------
187ClpNetworkMatrix::~ClpNetworkMatrix()
188{
189  delete matrix_;
190  delete[] lengths_;
191  delete[] indices_;
192}
193
194//----------------------------------------------------------------
195// Assignment operator
196//-------------------------------------------------------------------
197ClpNetworkMatrix &
198ClpNetworkMatrix::operator=(const ClpNetworkMatrix &rhs)
199{
200  if (this != &rhs) {
201    ClpMatrixBase::operator=(rhs);
202    delete matrix_;
203    delete[] lengths_;
204    delete[] indices_;
205    matrix_ = NULL;
206    lengths_ = NULL;
207    indices_ = NULL;
208    numberRows_ = rhs.numberRows_;
209    numberColumns_ = rhs.numberColumns_;
210    trueNetwork_ = rhs.trueNetwork_;
211    if (numberColumns_) {
212      indices_ = new int[2 * numberColumns_];
213      CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
214    }
215  }
216  return *this;
217}
218//-------------------------------------------------------------------
219// Clone
220//-------------------------------------------------------------------
221ClpMatrixBase *ClpNetworkMatrix::clone() const
222{
223  return new ClpNetworkMatrix(*this);
224}
225
226/* Returns a new matrix in reverse order without gaps */
227ClpMatrixBase *
228ClpNetworkMatrix::reverseOrderedCopy() const
229{
230  // count number in each row
231  CoinBigIndex *tempP = new CoinBigIndex[numberRows_];
232  CoinBigIndex *tempN = new CoinBigIndex[numberRows_];
233  memset(tempP, 0, numberRows_ * sizeof(CoinBigIndex));
234  memset(tempN, 0, numberRows_ * sizeof(CoinBigIndex));
235  CoinBigIndex j = 0;
236  int i;
237  for (i = 0; i < numberColumns_; i++, j += 2) {
238    int iRow = indices_[j];
239    tempN[iRow]++;
240    iRow = indices_[j + 1];
241    tempP[iRow]++;
242  }
243  int *newIndices = new int[2 * numberColumns_];
244  CoinBigIndex *newP = new CoinBigIndex[numberRows_ + 1];
245  CoinBigIndex *newN = new CoinBigIndex[numberRows_];
246  int iRow;
247  j = 0;
248  // do starts
249  for (iRow = 0; iRow < numberRows_; iRow++) {
250    newP[iRow] = j;
251    j += tempP[iRow];
252    tempP[iRow] = newP[iRow];
253    newN[iRow] = j;
254    j += tempN[iRow];
255    tempN[iRow] = newN[iRow];
256  }
257  newP[numberRows_] = j;
258  j = 0;
259  for (i = 0; i < numberColumns_; i++, j += 2) {
260    int iRow = indices_[j];
261    CoinBigIndex put = tempN[iRow];
262    newIndices[put++] = i;
263    tempN[iRow] = put;
264    iRow = indices_[j + 1];
265    put = tempP[iRow];
266    newIndices[put++] = i;
267    tempP[iRow] = put;
268  }
269  delete[] tempP;
270  delete[] tempN;
271  ClpPlusMinusOneMatrix *newCopy = new ClpPlusMinusOneMatrix();
272  newCopy->passInCopy(numberRows_, numberColumns_,
273    false, newIndices, newP, newN);
274  return newCopy;
275}
276//unscaled versions
277void ClpNetworkMatrix::times(double scalar,
278  const double *x, double *y) const
279{
280  int iColumn;
281  CoinBigIndex j = 0;
282  if (trueNetwork_) {
283    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
284      double value = scalar * x[iColumn];
285      if (value) {
286        int iRowM = indices_[j];
287        int iRowP = indices_[j + 1];
288        y[iRowM] -= value;
289        y[iRowP] += value;
290      }
291    }
292  } else {
293    // skip negative rows
294    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
295      double value = scalar * x[iColumn];
296      if (value) {
297        int iRowM = indices_[j];
298        int iRowP = indices_[j + 1];
299        if (iRowM >= 0)
300          y[iRowM] -= value;
301        if (iRowP >= 0)
302          y[iRowP] += value;
303      }
304    }
305  }
306}
307void ClpNetworkMatrix::transposeTimes(double scalar,
308  const double *x, double *y) const
309{
310  int iColumn;
311  CoinBigIndex j = 0;
312  if (trueNetwork_) {
313    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
314      double value = y[iColumn];
315      int iRowM = indices_[j];
316      int iRowP = indices_[j + 1];
317      value -= scalar * x[iRowM];
318      value += scalar * x[iRowP];
319      y[iColumn] = value;
320    }
321  } else {
322    // skip negative rows
323    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
324      double value = y[iColumn];
325      int iRowM = indices_[j];
326      int iRowP = indices_[j + 1];
327      if (iRowM >= 0)
328        value -= scalar * x[iRowM];
329      if (iRowP >= 0)
330        value += scalar * x[iRowP];
331      y[iColumn] = value;
332    }
333  }
334}
335void ClpNetworkMatrix::times(double scalar,
336  const double *x, double *y,
337  const double * /*rowScale*/,
338  const double * /*columnScale*/) const
339{
340  // we know it is not scaled
341  times(scalar, x, y);
342}
343void ClpNetworkMatrix::transposeTimes(double scalar,
344  const double *x, double *y,
345  const double * /*rowScale*/,
346  const double * /*columnScale*/,
347  double * /*spare*/) const
348{
349  // we know it is not scaled
350  transposeTimes(scalar, x, y);
351}
352/* Return <code>x * A + y</code> in <code>z</code>.
353        Squashes small elements and knows about ClpSimplex */
354void ClpNetworkMatrix::transposeTimes(const ClpSimplex *model, double scalar,
355  const CoinIndexedVector *rowArray,
356  CoinIndexedVector *y,
357  CoinIndexedVector *columnArray) const
358{
359  // we know it is not scaled
360  columnArray->clear();
361  double *pi = rowArray->denseVector();
362  int numberNonZero = 0;
363  int *index = columnArray->getIndices();
364  double *array = columnArray->denseVector();
365  int numberInRowArray = rowArray->getNumElements();
366  // maybe I need one in OsiSimplex
367  double zeroTolerance = model->zeroTolerance();
368  int numberRows = model->numberRows();
369#ifndef NO_RTTI
370  ClpPlusMinusOneMatrix *rowCopy = dynamic_cast< ClpPlusMinusOneMatrix * >(model->rowCopy());
371#else
372  ClpPlusMinusOneMatrix *rowCopy = static_cast< ClpPlusMinusOneMatrix * >(model->rowCopy());
373#endif
374  bool packed = rowArray->packedMode();
375  double factor = 0.3;
376  // We may not want to do by row if there may be cache problems
377  int numberColumns = model->numberColumns();
378  // It would be nice to find L2 cache size - for moment 512K
379  // Be slightly optimistic
380  if (numberColumns * sizeof(double) > 1000000) {
381    if (numberRows * 10 < numberColumns)
382      factor = 0.1;
383    else if (numberRows * 4 < numberColumns)
384      factor = 0.15;
385    else if (numberRows * 2 < numberColumns)
386      factor = 0.2;
387    //if (model->numberIterations()%50==0)
388    //printf("%d nonzero\n",numberInRowArray);
389  }
390  if (numberInRowArray > factor * numberRows || !rowCopy) {
391    // do by column
392    int iColumn;
393    assert(!y->getNumElements());
394    CoinBigIndex j = 0;
395    if (packed) {
396      // need to expand pi into y
397      assert(y->capacity() >= numberRows);
398      double *piOld = pi;
399      pi = y->denseVector();
400      const int *whichRow = rowArray->getIndices();
401      int i;
402      // modify pi so can collapse to one loop
403      for (i = 0; i < numberInRowArray; i++) {
404        int iRow = whichRow[i];
405        pi[iRow] = scalar * piOld[i];
406      }
407      if (trueNetwork_) {
408        for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
409          double value = 0.0;
410          int iRowM = indices_[j];
411          int iRowP = indices_[j + 1];
412          value -= pi[iRowM];
413          value += pi[iRowP];
414          if (fabs(value) > zeroTolerance) {
415            array[numberNonZero] = value;
416            index[numberNonZero++] = iColumn;
417          }
418        }
419      } else {
420        // skip negative rows
421        for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
422          double value = 0.0;
423          int iRowM = indices_[j];
424          int iRowP = indices_[j + 1];
425          if (iRowM >= 0)
426            value -= pi[iRowM];
427          if (iRowP >= 0)
428            value += pi[iRowP];
429          if (fabs(value) > zeroTolerance) {
430            array[numberNonZero] = value;
431            index[numberNonZero++] = iColumn;
432          }
433        }
434      }
435      for (i = 0; i < numberInRowArray; i++) {
436        int iRow = whichRow[i];
437        pi[iRow] = 0.0;
438      }
439    } else {
440      if (trueNetwork_) {
441        for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
442          double value = 0.0;
443          int iRowM = indices_[j];
444          int iRowP = indices_[j + 1];
445          value -= scalar * pi[iRowM];
446          value += scalar * pi[iRowP];
447          if (fabs(value) > zeroTolerance) {
448            index[numberNonZero++] = iColumn;
449            array[iColumn] = value;
450          }
451        }
452      } else {
453        // skip negative rows
454        for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
455          double value = 0.0;
456          int iRowM = indices_[j];
457          int iRowP = indices_[j + 1];
458          if (iRowM >= 0)
459            value -= scalar * pi[iRowM];
460          if (iRowP >= 0)
461            value += scalar * pi[iRowP];
462          if (fabs(value) > zeroTolerance) {
463            index[numberNonZero++] = iColumn;
464            array[iColumn] = value;
465          }
466        }
467      }
468    }
469    columnArray->setNumElements(numberNonZero);
470  } else {
471    // do by row
472    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
473  }
474}
475/* Return <code>x *A in <code>z</code> but
476   just for indices in y. */
477void ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/,
478  const CoinIndexedVector *rowArray,
479  const CoinIndexedVector *y,
480  CoinIndexedVector *columnArray) const
481{
482  columnArray->clear();
483  double *pi = rowArray->denseVector();
484  double *array = columnArray->denseVector();
485  int jColumn;
486  int numberToDo = y->getNumElements();
487  const int *which = y->getIndices();
488  assert(!rowArray->packedMode());
489  columnArray->setPacked();
490  if (trueNetwork_) {
491    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
492      int iColumn = which[jColumn];
493      double value = 0.0;
494      CoinBigIndex j = iColumn << 1;
495      int iRowM = indices_[j];
496      int iRowP = indices_[j + 1];
497      value -= pi[iRowM];
498      value += pi[iRowP];
499      array[jColumn] = value;
500    }
501  } else {
502    // skip negative rows
503    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
504      int iColumn = which[jColumn];
505      double value = 0.0;
506      CoinBigIndex j = iColumn << 1;
507      int iRowM = indices_[j];
508      int iRowP = indices_[j + 1];
509      if (iRowM >= 0)
510        value -= pi[iRowM];
511      if (iRowP >= 0)
512        value += pi[iRowP];
513      array[jColumn] = value;
514    }
515  }
516}
517/// returns number of elements in column part of basis,
518int ClpNetworkMatrix::countBasis(const int *whichColumn,
519  int &numberColumnBasic)
520{
521  int i;
522  CoinBigIndex numberElements = 0;
523  if (trueNetwork_) {
524    numberElements = 2 * numberColumnBasic;
525  } else {
526    for (i = 0; i < numberColumnBasic; i++) {
527      int iColumn = whichColumn[i];
528      CoinBigIndex j = iColumn << 1;
529      int iRowM = indices_[j];
530      int iRowP = indices_[j + 1];
531      if (iRowM >= 0)
532        numberElements++;
533      if (iRowP >= 0)
534        numberElements++;
535    }
536  }
537#if COIN_BIG_INDEX
538  if (numberElements > COIN_INT_MAX) {
539    printf("Factorization too large\n");
540    abort();
541  }
542#endif
543  return static_cast< int >(numberElements);
544}
545void ClpNetworkMatrix::fillBasis(ClpSimplex * /*model*/,
546  const int *whichColumn,
547  int &numberColumnBasic,
548  int *indexRowU, int *start,
549  int *rowCount, int *columnCount,
550  CoinFactorizationDouble *elementU)
551{
552  int i;
553  CoinBigIndex numberElements = start[0];
554  if (trueNetwork_) {
555    for (i = 0; i < numberColumnBasic; i++) {
556      int iColumn = whichColumn[i];
557      CoinBigIndex j = iColumn << 1;
558      int iRowM = indices_[j];
559      int iRowP = indices_[j + 1];
560      indexRowU[numberElements] = iRowM;
561      rowCount[iRowM]++;
562      elementU[numberElements] = -1.0;
563      indexRowU[numberElements + 1] = iRowP;
564      rowCount[iRowP]++;
565      elementU[numberElements + 1] = 1.0;
566      numberElements += 2;
567      start[i + 1] = static_cast< int >(numberElements);
568      columnCount[i] = 2;
569    }
570  } else {
571    for (i = 0; i < numberColumnBasic; i++) {
572      int iColumn = whichColumn[i];
573      CoinBigIndex j = iColumn << 1;
574      int iRowM = indices_[j];
575      int iRowP = indices_[j + 1];
576      if (iRowM >= 0) {
577        indexRowU[numberElements] = iRowM;
578        rowCount[iRowM]++;
579        elementU[numberElements++] = -1.0;
580      }
581      if (iRowP >= 0) {
582        indexRowU[numberElements] = iRowP;
583        rowCount[iRowP]++;
584        elementU[numberElements++] = 1.0;
585      }
586      start[i + 1] = static_cast< int >(numberElements);
587      columnCount[i] = static_cast< int >(numberElements) - start[i];
588    }
589  }
590#if COIN_BIG_INDEX
591  if (numberElements > COIN_INT_MAX)
592    abort();
593#endif
594}
595/* Unpacks a column into an CoinIndexedvector
596 */
597void ClpNetworkMatrix::unpack(const ClpSimplex * /*model*/, CoinIndexedVector *rowArray,
598  int iColumn) const
599{
600  CoinBigIndex j = iColumn << 1;
601  int iRowM = indices_[j];
602  int iRowP = indices_[j + 1];
603  if (iRowM >= 0)
604    rowArray->add(iRowM, -1.0);
605  if (iRowP >= 0)
606    rowArray->add(iRowP, 1.0);
607}
608/* Unpacks a column into an CoinIndexedvector
609** in packed foramt
610Note that model is NOT const.  Bounds and objective could
611be modified if doing column generation (just for this variable) */
612void ClpNetworkMatrix::unpackPacked(ClpSimplex * /*model*/,
613  CoinIndexedVector *rowArray,
614  int iColumn) const
615{
616  int *index = rowArray->getIndices();
617  double *array = rowArray->denseVector();
618  int number = 0;
619  CoinBigIndex j = iColumn << 1;
620  int iRowM = indices_[j];
621  int iRowP = indices_[j + 1];
622  if (iRowM >= 0) {
623    array[number] = -1.0;
624    index[number++] = iRowM;
625  }
626  if (iRowP >= 0) {
627    array[number] = 1.0;
628    index[number++] = iRowP;
629  }
630  rowArray->setNumElements(number);
631  rowArray->setPackedMode(true);
632}
633/* Adds multiple of a column into an CoinIndexedvector
634      You can use quickAdd to add to vector */
635void ClpNetworkMatrix::add(const ClpSimplex * /*model*/, CoinIndexedVector *rowArray,
636  int iColumn, double multiplier) const
637{
638  CoinBigIndex j = iColumn << 1;
639  int iRowM = indices_[j];
640  int iRowP = indices_[j + 1];
641  if (iRowM >= 0)
642    rowArray->quickAdd(iRowM, -multiplier);
643  if (iRowP >= 0)
644    rowArray->quickAdd(iRowP, multiplier);
645}
646/* Adds multiple of a column into an array */
647void ClpNetworkMatrix::add(const ClpSimplex * /*model*/, double *array,
648  int iColumn, double multiplier) const
649{
650  CoinBigIndex j = iColumn << 1;
651  int iRowM = indices_[j];
652  int iRowP = indices_[j + 1];
653  if (iRowM >= 0)
654    array[iRowM] -= multiplier;
655  if (iRowP >= 0)
656    array[iRowP] += multiplier;
657}
658
659// Return a complete CoinPackedMatrix
660CoinPackedMatrix *
661ClpNetworkMatrix::getPackedMatrix() const
662{
663  if (!matrix_) {
664    assert(trueNetwork_); // fix later
665    int numberElements = 2 * numberColumns_;
666    double *elements = new double[numberElements];
667    CoinBigIndex i;
668    for (i = 0; i < 2 * numberColumns_; i += 2) {
669      elements[i] = -1.0;
670      elements[i + 1] = 1.0;
671    }
672    CoinBigIndex *starts = new CoinBigIndex[numberColumns_ + 1];
673    for (i = 0; i < numberColumns_ + 1; i++) {
674      starts[i] = 2 * i;
675    }
676    // use assignMatrix to save space
677    delete[] lengths_;
678    lengths_ = NULL;
679    matrix_ = new CoinPackedMatrix();
680    int *indices = CoinCopyOfArray(indices_, 2 * numberColumns_);
681    matrix_->assignMatrix(true, numberRows_, numberColumns_,
682      getNumElements(),
683      elements, indices,
684      starts, lengths_);
685    assert(!elements);
686    assert(!starts);
687    assert(!indices);
688    assert(!lengths_);
689  }
690  return matrix_;
691}
692/* A vector containing the elements in the packed matrix. Note that there
693   might be gaps in this list, entries that do not belong to any
694   major-dimension vector. To get the actual elements one should look at
695   this vector together with vectorStarts and vectorLengths. */
696const double *
697ClpNetworkMatrix::getElements() const
698{
699  if (!matrix_)
700    getPackedMatrix();
701  return matrix_->getElements();
702}
703
704const CoinBigIndex *
705ClpNetworkMatrix::getVectorStarts() const
706{
707  if (!matrix_)
708    getPackedMatrix();
709  return matrix_->getVectorStarts();
710}
711/* The lengths of the major-dimension vectors. */
712const int *
713ClpNetworkMatrix::getVectorLengths() const
714{
715  assert(trueNetwork_); // fix later
716  if (!lengths_) {
717    lengths_ = new int[numberColumns_];
718    int i;
719    for (i = 0; i < numberColumns_; i++) {
720      lengths_[i] = 2;
721    }
722  }
723  return lengths_;
724}
725/* Delete the columns whose indices are listed in <code>indDel</code>. */
726void ClpNetworkMatrix::deleteCols(const int numDel, const int *indDel)
727{
728  assert(trueNetwork_);
729  int iColumn;
730  int numberBad = 0;
731  // Use array to make sure we can have duplicates
732  char *which = new char[numberColumns_];
733  memset(which, 0, numberColumns_);
734  int nDuplicate = 0;
735  for (iColumn = 0; iColumn < numDel; iColumn++) {
736    int jColumn = indDel[iColumn];
737    if (jColumn < 0 || jColumn >= numberColumns_) {
738      numberBad++;
739    } else {
740      if (which[jColumn])
741        nDuplicate++;
742      else
743        which[jColumn] = 1;
744    }
745  }
746  if (numberBad)
747    throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
748  int newNumber = numberColumns_ - numDel + nDuplicate;
749  // Get rid of temporary arrays
750  delete[] lengths_;
751  lengths_ = NULL;
752  delete matrix_;
753  matrix_ = NULL;
754  int newSize = 2 * newNumber;
755  int *newIndices = new int[newSize];
756  newSize = 0;
757  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
758    if (!which[iColumn]) {
759      CoinBigIndex start;
760      CoinBigIndex i;
761      start = 2 * iColumn;
762      for (i = start; i < start + 2; i++)
763        newIndices[newSize++] = indices_[i];
764    }
765  }
766  delete[] which;
767  delete[] indices_;
768  indices_ = newIndices;
769  numberColumns_ = newNumber;
770}
771/* Delete the rows whose indices are listed in <code>indDel</code>. */
772void ClpNetworkMatrix::deleteRows(const int numDel, const int *indDel)
773{
774  int iRow;
775  int numberBad = 0;
776  // Use array to make sure we can have duplicates
777  int *which = new int[numberRows_];
778  memset(which, 0, numberRows_ * sizeof(int));
779  for (iRow = 0; iRow < numDel; iRow++) {
780    int jRow = indDel[iRow];
781    if (jRow < 0 || jRow >= numberRows_) {
782      numberBad++;
783    } else {
784      which[jRow] = 1;
785    }
786  }
787  if (numberBad)
788    throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
789  // Only valid of all columns have 0 entries
790  int iColumn;
791  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
792    CoinBigIndex start;
793    CoinBigIndex i;
794    start = 2 * iColumn;
795    for (i = start; i < start + 2; i++) {
796      int iRow = indices_[i];
797      if (which[iRow])
798        numberBad++;
799    }
800  }
801  if (numberBad)
802    throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
803  int newNumber = 0;
804  for (iRow = 0; iRow < numberRows_; iRow++) {
805    if (!which[iRow])
806      which[iRow] = newNumber++;
807    else
808      which[iRow] = -1;
809  }
810  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
811    CoinBigIndex start;
812    CoinBigIndex i;
813    start = 2 * iColumn;
814    for (i = start; i < start + 2; i++) {
815      int iRow = indices_[i];
816      indices_[i] = which[iRow];
817    }
818  }
819  delete[] which;
820  numberRows_ = newNumber;
821}
822/* Given positive integer weights for each row fills in sum of weights
823   for each column (and slack).
824   Returns weights vector
825*/
826CoinBigIndex *
827ClpNetworkMatrix::dubiousWeights(const ClpSimplex *model, int *inputWeights) const
828{
829  int numberRows = model->numberRows();
830  int numberColumns = model->numberColumns();
831  int number = numberRows + numberColumns;
832  CoinBigIndex *weights = new CoinBigIndex[number];
833  int i;
834  for (i = 0; i < numberColumns; i++) {
835    CoinBigIndex j = i << 1;
836    CoinBigIndex count = 0;
837    int iRowM = indices_[j];
838    int iRowP = indices_[j + 1];
839    if (iRowM >= 0) {
840      count += inputWeights[iRowM];
841    }
842    if (iRowP >= 0) {
843      count += inputWeights[iRowP];
844    }
845    weights[i] = count;
846  }
847  for (i = 0; i < numberRows; i++) {
848    weights[i + numberColumns] = inputWeights[i];
849  }
850  return weights;
851}
852/* Returns largest and smallest elements of both signs.
853   Largest refers to largest absolute value.
854*/
855void ClpNetworkMatrix::rangeOfElements(double &smallestNegative, double &largestNegative,
856  double &smallestPositive, double &largestPositive)
857{
858  smallestNegative = -1.0;
859  largestNegative = -1.0;
860  smallestPositive = 1.0;
861  largestPositive = 1.0;
862}
863// Says whether it can do partial pricing
864bool ClpNetworkMatrix::canDoPartialPricing() const
865{
866  return true;
867}
868// Partial pricing
869void ClpNetworkMatrix::partialPricing(ClpSimplex *model, double startFraction, double endFraction,
870  int &bestSequence, int &numberWanted)
871{
872  numberWanted = currentWanted_;
873  int j;
874  int start = static_cast< int >(startFraction * numberColumns_);
875  int end = CoinMin(static_cast< int >(endFraction * numberColumns_ + 1), numberColumns_);
876  double tolerance = model->currentDualTolerance();
877  double *reducedCost = model->djRegion();
878  const double *duals = model->dualRowSolution();
879  const double *cost = model->costRegion();
880  double bestDj;
881  if (bestSequence >= 0)
882    bestDj = fabs(reducedCost[bestSequence]);
883  else
884    bestDj = tolerance;
885  int sequenceOut = model->sequenceOut();
886  int saveSequence = bestSequence;
887  if (!trueNetwork_) {
888    // Not true network
889    int iSequence;
890    for (iSequence = start; iSequence < end; iSequence++) {
891      if (iSequence != sequenceOut) {
892        double value;
893        int iRowM, iRowP;
894        ClpSimplex::Status status = model->getStatus(iSequence);
895
896        switch (status) {
897
898        case ClpSimplex::basic:
899        case ClpSimplex::isFixed:
900          break;
901        case ClpSimplex::isFree:
902        case ClpSimplex::superBasic:
903          value = cost[iSequence];
904          j = iSequence << 1;
905          // skip negative rows
906          iRowM = indices_[j];
907          iRowP = indices_[j + 1];
908          if (iRowM >= 0)
909            value += duals[iRowM];
910          if (iRowP >= 0)
911            value -= duals[iRowP];
912          value = fabs(value);
913          if (value > FREE_ACCEPT * tolerance) {
914            numberWanted--;
915            // we are going to bias towards free (but only if reasonable)
916            value *= FREE_BIAS;
917            if (value > bestDj) {
918              // check flagged variable and correct dj
919              if (!model->flagged(iSequence)) {
920                bestDj = value;
921                bestSequence = iSequence;
922              } else {
923                // just to make sure we don't exit before got something
924                numberWanted++;
925              }
926            }
927          }
928          break;
929        case ClpSimplex::atUpperBound:
930          value = cost[iSequence];
931          j = iSequence << 1;
932          // skip negative rows
933          iRowM = indices_[j];
934          iRowP = indices_[j + 1];
935          if (iRowM >= 0)
936            value += duals[iRowM];
937          if (iRowP >= 0)
938            value -= duals[iRowP];
939          if (value > tolerance) {
940            numberWanted--;
941            if (value > bestDj) {
942              // check flagged variable and correct dj
943              if (!model->flagged(iSequence)) {
944                bestDj = value;
945                bestSequence = iSequence;
946              } else {
947                // just to make sure we don't exit before got something
948                numberWanted++;
949              }
950            }
951          }
952          break;
953        case ClpSimplex::atLowerBound:
954          value = cost[iSequence];
955          j = iSequence << 1;
956          // skip negative rows
957          iRowM = indices_[j];
958          iRowP = indices_[j + 1];
959          if (iRowM >= 0)
960            value += duals[iRowM];
961          if (iRowP >= 0)
962            value -= duals[iRowP];
963          value = -value;
964          if (value > tolerance) {
965            numberWanted--;
966            if (value > bestDj) {
967              // check flagged variable and correct dj
968              if (!model->flagged(iSequence)) {
969                bestDj = value;
970                bestSequence = iSequence;
971              } else {
972                // just to make sure we don't exit before got something
973                numberWanted++;
974              }
975            }
976          }
977          break;
978        }
979      }
980      if (!numberWanted)
981        break;
982    }
983    if (bestSequence != saveSequence) {
984      // recompute dj
985      double value = cost[bestSequence];
986      j = bestSequence << 1;
987      // skip negative rows
988      int iRowM = indices_[j];
989      int iRowP = indices_[j + 1];
990      if (iRowM >= 0)
991        value += duals[iRowM];
992      if (iRowP >= 0)
993        value -= duals[iRowP];
994      reducedCost[bestSequence] = value;
995      savedBestSequence_ = bestSequence;
996      savedBestDj_ = reducedCost[savedBestSequence_];
997    }
998  } else {
999    // true network
1000    int iSequence;
1001    for (iSequence = start; iSequence < end; iSequence++) {
1002      if (iSequence != sequenceOut) {
1003        double value;
1004        int iRowM, iRowP;
1005        ClpSimplex::Status status = model->getStatus(iSequence);
1006
1007        switch (status) {
1008
1009        case ClpSimplex::basic:
1010        case ClpSimplex::isFixed:
1011          break;
1012        case ClpSimplex::isFree:
1013        case ClpSimplex::superBasic:
1014          value = cost[iSequence];
1015          j = iSequence << 1;
1016          iRowM = indices_[j];
1017          iRowP = indices_[j + 1];
1018          value += duals[iRowM];
1019          value -= duals[iRowP];
1020          value = fabs(value);
1021          if (value > FREE_ACCEPT * tolerance) {
1022            numberWanted--;
1023            // we are going to bias towards free (but only if reasonable)
1024            value *= FREE_BIAS;
1025            if (value > bestDj) {
1026              // check flagged variable and correct dj
1027              if (!model->flagged(iSequence)) {
1028                bestDj = value;
1029                bestSequence = iSequence;
1030              } else {
1031                // just to make sure we don't exit before got something
1032                numberWanted++;
1033              }
1034            }
1035          }
1036          break;
1037        case ClpSimplex::atUpperBound:
1038          value = cost[iSequence];
1039          j = iSequence << 1;
1040          iRowM = indices_[j];
1041          iRowP = indices_[j + 1];
1042          value += duals[iRowM];
1043          value -= duals[iRowP];
1044          if (value > tolerance) {
1045            numberWanted--;
1046            if (value > bestDj) {
1047              // check flagged variable and correct dj
1048              if (!model->flagged(iSequence)) {
1049                bestDj = value;
1050                bestSequence = iSequence;
1051              } else {
1052                // just to make sure we don't exit before got something
1053                numberWanted++;
1054              }
1055            }
1056          }
1057          break;
1058        case ClpSimplex::atLowerBound:
1059          value = cost[iSequence];
1060          j = iSequence << 1;
1061          iRowM = indices_[j];
1062          iRowP = indices_[j + 1];
1063          value += duals[iRowM];
1064          value -= duals[iRowP];
1065          value = -value;
1066          if (value > tolerance) {
1067            numberWanted--;
1068            if (value > bestDj) {
1069              // check flagged variable and correct dj
1070              if (!model->flagged(iSequence)) {
1071                bestDj = value;
1072                bestSequence = iSequence;
1073              } else {
1074                // just to make sure we don't exit before got something
1075                numberWanted++;
1076              }
1077            }
1078          }
1079          break;
1080        }
1081      }
1082      if (!numberWanted)
1083        break;
1084    }
1085    if (bestSequence != saveSequence) {
1086      // recompute dj
1087      double value = cost[bestSequence];
1088      j = bestSequence << 1;
1089      int iRowM = indices_[j];
1090      int iRowP = indices_[j + 1];
1091      value += duals[iRowM];
1092      value -= duals[iRowP];
1093      reducedCost[bestSequence] = value;
1094      savedBestSequence_ = bestSequence;
1095      savedBestDj_ = reducedCost[savedBestSequence_];
1096    }
1097  }
1098  currentWanted_ = numberWanted;
1099}
1100// Allow any parts of a created CoinMatrix to be deleted
1101void ClpNetworkMatrix::releasePackedMatrix() const
1102{
1103  delete matrix_;
1104  delete[] lengths_;
1105  matrix_ = NULL;
1106  lengths_ = NULL;
1107}
1108// Append Columns
1109void ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase *const *columns)
1110{
1111  int iColumn;
1112  int numberBad = 0;
1113  for (iColumn = 0; iColumn < number; iColumn++) {
1114    int n = columns[iColumn]->getNumElements();
1115    const double *element = columns[iColumn]->getElements();
1116    if (n != 2)
1117      numberBad++;
1118    if (fabs(element[0]) != 1.0 || fabs(element[1]) != 1.0)
1119      numberBad++;
1120    else if (element[0] * element[1] != -1.0)
1121      numberBad++;
1122  }
1123  if (numberBad)
1124    throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
1125  // Get rid of temporary arrays
1126  delete[] lengths_;
1127  lengths_ = NULL;
1128  delete matrix_;
1129  matrix_ = NULL;
1130  CoinBigIndex size = 2 * number;
1131  int *temp2 = new int[numberColumns_ * 2 + size];
1132  CoinMemcpyN(indices_, numberColumns_ * 2, temp2);
1133  delete[] indices_;
1134  indices_ = temp2;
1135  // now add
1136  size = 2 * numberColumns_;
1137  for (iColumn = 0; iColumn < number; iColumn++) {
1138    const int *row = columns[iColumn]->getIndices();
1139    const double *element = columns[iColumn]->getElements();
1140    if (element[0] == -1.0) {
1141      indices_[size++] = row[0];
1142      indices_[size++] = row[1];
1143    } else {
1144      indices_[size++] = row[1];
1145      indices_[size++] = row[0];
1146    }
1147  }
1148
1149  numberColumns_ += number;
1150}
1151// Append Rows
1152void ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase *const *rows)
1153{
1154  // must be zero arrays
1155  int numberBad = 0;
1156  int iRow;
1157  for (iRow = 0; iRow < number; iRow++) {
1158    numberBad += rows[iRow]->getNumElements();
1159  }
1160  if (numberBad)
1161    throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
1162  numberRows_ += number;
1163}
1164#ifndef SLIM_CLP
1165/* Append a set of rows/columns to the end of the matrix. Returns number of errors
1166   i.e. if any of the new rows/columns contain an index that's larger than the
1167   number of columns-1/rows-1 (if numberOther>0) or duplicates
1168   If 0 then rows, 1 if columns */
1169int ClpNetworkMatrix::appendMatrix(int number, int type,
1170  const CoinBigIndex *starts, const int *index,
1171  const double *element, int /*numberOther*/)
1172{
1173  int numberErrors = 0;
1174  // make into CoinPackedVector
1175  CoinPackedVectorBase **vectors = new CoinPackedVectorBase *[number];
1176  int iVector;
1177  for (iVector = 0; iVector < number; iVector++) {
1178    CoinBigIndex iStart = starts[iVector];
1179    vectors[iVector] = new CoinPackedVector(static_cast< int >(starts[iVector + 1] - iStart),
1180      index + iStart, element + iStart);
1181  }
1182  if (type == 0) {
1183    // rows
1184    appendRows(number, vectors);
1185  } else {
1186    // columns
1187    appendCols(number, vectors);
1188  }
1189  for (iVector = 0; iVector < number; iVector++)
1190    delete vectors[iVector];
1191  delete[] vectors;
1192  return numberErrors;
1193}
1194#endif
1195/* Subset clone (without gaps).  Duplicates are allowed
1196   and order is as given */
1197ClpMatrixBase *
1198ClpNetworkMatrix::subsetClone(int numberRows, const int *whichRows,
1199  int numberColumns,
1200  const int *whichColumns) const
1201{
1202  return new ClpNetworkMatrix(*this, numberRows, whichRows,
1203    numberColumns, whichColumns);
1204}
1205/* Subset constructor (without gaps).  Duplicates are allowed
1206   and order is as given */
1207ClpNetworkMatrix::ClpNetworkMatrix(
1208  const ClpNetworkMatrix &rhs,
1209  int numberRows, const int *whichRow,
1210  int numberColumns, const int *whichColumn)
1211  : ClpMatrixBase(rhs)
1212{
1213  setType(11);
1214  matrix_ = NULL;
1215  lengths_ = NULL;
1216  indices_ = new int[2 * numberColumns];
1217  ;
1218  numberRows_ = numberRows;
1219  numberColumns_ = numberColumns;
1220  trueNetwork_ = true;
1221  int iColumn;
1222  int numberBad = 0;
1223  int *which = new int[rhs.numberRows_];
1224  int iRow;
1225  for (iRow = 0; iRow < rhs.numberRows_; iRow++)
1226    which[iRow] = -1;
1227  int n = 0;
1228  for (iRow = 0; iRow < numberRows; iRow++) {
1229    int jRow = whichRow[iRow];
1230    assert(jRow >= 0 && jRow < rhs.numberRows_);
1231    which[jRow] = n++;
1232  }
1233  for (iColumn = 0; iColumn < numberColumns; iColumn++) {
1234    CoinBigIndex start;
1235    CoinBigIndex i;
1236    start = 2 * iColumn;
1237    CoinBigIndex offset = 2 * whichColumn[iColumn] - start;
1238    for (i = start; i < start + 2; i++) {
1239      int iRow = rhs.indices_[i + offset];
1240      iRow = which[iRow];
1241      if (iRow < 0)
1242        numberBad++;
1243      else
1244        indices_[i] = iRow;
1245    }
1246  }
1247  if (numberBad)
1248    throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
1249}
1250
1251/* vi: softtabstop=2 shiftwidth=2 expandtab tabstop=2
1252*/
Note: See TracBrowser for help on using the repository browser.