source: trunk/Clp/src/ClpNetworkMatrix.cpp @ 1370

Last change on this file since 1370 was 1370, checked in by forrest, 11 years ago

add ids

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