source: branches/devel/Clp/src/ClpNetworkMatrix.cpp @ 929

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

fix bug

  • Property svn:eol-style set to native
  • Property svn:keywords set to Author Date Id Revision
File size: 26.7 KB
Line 
1// Copyright (C) 2003, International Business Machines
2// Corporation and others.  All Rights Reserved.
3
4
5#include <cstdio>
6
7#include "CoinPragma.hpp"
8#include "CoinIndexedVector.hpp"
9#include "CoinHelperFunctions.hpp"
10
11#include "ClpSimplex.hpp"
12#include "ClpFactorization.hpp"
13// at end to get min/max!
14#include "ClpNetworkMatrix.hpp"
15#include "ClpPlusMinusOneMatrix.hpp"
16#include "ClpMessage.hpp"
17#include <iostream>
18
19//#############################################################################
20// Constructors / Destructor / Assignment
21//#############################################################################
22
23//-------------------------------------------------------------------
24// Default Constructor
25//-------------------------------------------------------------------
26ClpNetworkMatrix::ClpNetworkMatrix () 
27  : ClpMatrixBase()
28{
29  setType(11);
30  elements_ = NULL;
31  starts_ = NULL;
32  lengths_=NULL;
33  indices_=NULL;
34  numberRows_=0;
35  numberColumns_=0;
36  trueNetwork_=false;
37}
38
39/* Constructor from two arrays */
40ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int * head,
41                                   const int * tail)
42  : ClpMatrixBase()
43{
44  setType(11);
45  elements_ = NULL;
46  starts_ = NULL;
47  lengths_=NULL;
48  indices_=new int[2*numberColumns];;
49  numberRows_=-1;
50  numberColumns_=numberColumns;
51  trueNetwork_=true;
52  int iColumn;
53  CoinBigIndex j=0;
54  for (iColumn=0;iColumn<numberColumns_;iColumn++, j+=2) {
55    int iRow = head[iColumn];
56    numberRows_ = CoinMax(numberRows_,iRow);
57    indices_[j]=iRow;
58    iRow = tail[iColumn];
59    numberRows_ = CoinMax(numberRows_,iRow);
60    indices_[j+1]=iRow;
61  }
62  numberRows_++;
63}
64//-------------------------------------------------------------------
65// Copy constructor
66//-------------------------------------------------------------------
67ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix & rhs) 
68: ClpMatrixBase(rhs)
69{ 
70  elements_ = NULL;
71  starts_ = 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    memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
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  elements_ = NULL;
94  starts_ = 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 [] elements_;
190  delete [] starts_;
191  delete [] lengths_;
192  delete [] indices_;
193}
194
195//----------------------------------------------------------------
196// Assignment operator
197//-------------------------------------------------------------------
198ClpNetworkMatrix &
199ClpNetworkMatrix::operator=(const ClpNetworkMatrix& rhs)
200{
201  if (this != &rhs) {
202    ClpMatrixBase::operator=(rhs);
203    delete [] elements_;
204    delete [] starts_;
205    delete [] lengths_;
206    delete [] indices_;
207    elements_ = NULL;
208    starts_ = NULL;
209    lengths_=NULL;
210    indices_=NULL;
211    numberRows_=rhs.numberRows_;
212    numberColumns_=rhs.numberColumns_;
213    trueNetwork_=rhs.trueNetwork_;
214    if (numberColumns_) {
215      indices_ = new int [ 2*numberColumns_];
216      memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
217    }
218  }
219  return *this;
220}
221//-------------------------------------------------------------------
222// Clone
223//-------------------------------------------------------------------
224ClpMatrixBase * ClpNetworkMatrix::clone() const
225{
226  return new ClpNetworkMatrix(*this);
227}
228
229/* Returns a new matrix in reverse order without gaps */
230ClpMatrixBase * 
231ClpNetworkMatrix::reverseOrderedCopy() const
232{
233  // count number in each row
234  CoinBigIndex * tempP = new CoinBigIndex [numberRows_];
235  CoinBigIndex * tempN = new CoinBigIndex [numberRows_];
236  memset(tempP,0,numberRows_*sizeof(CoinBigIndex));
237  memset(tempN,0,numberRows_*sizeof(CoinBigIndex));
238  CoinBigIndex j=0;
239  int i;
240  for (i=0;i<numberColumns_;i++,j+=2) {
241    int iRow = indices_[j];
242    tempN[iRow]++;
243    iRow = indices_[j+1];
244    tempP[iRow]++;
245  }
246  int * newIndices = new int [2*numberColumns_];
247  CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];
248  CoinBigIndex * newN = new CoinBigIndex[numberRows_];
249  int iRow;
250  j=0;
251  // do starts
252  for (iRow=0;iRow<numberRows_;iRow++) {
253    newP[iRow]=j;
254    j += tempP[iRow];
255    tempP[iRow]=newP[iRow];
256    newN[iRow] = j;
257    j += tempN[iRow];
258    tempN[iRow]=newN[iRow];
259  }
260  newP[numberRows_]=j;
261  j=0;
262  for (i=0;i<numberColumns_;i++,j+=2) {
263    int iRow = indices_[j];
264    CoinBigIndex put = tempN[iRow];
265    newIndices[put++] = i;
266    tempN[iRow] = put;
267    iRow = indices_[j+1];
268    put = tempP[iRow];
269    newIndices[put++] = i;
270    tempP[iRow] = put;
271  }
272  delete [] tempP;
273  delete [] tempN;
274  ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
275  newCopy->passInCopy(numberRows_, numberColumns_,
276                      false,  newIndices, newP, newN);
277  return newCopy;
278}
279//unscaled versions
280void 
281ClpNetworkMatrix::times(double scalar,
282                   const double * x, double * y) const
283{
284  int iColumn;
285  CoinBigIndex j=0;
286  if (trueNetwork_) {
287    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
288      double value = scalar*x[iColumn];
289      if (value) {
290        int iRowM = indices_[j];
291        int iRowP = indices_[j+1];
292        y[iRowM] -= value;
293        y[iRowP] += value;
294      }
295    }
296  } else {
297    // skip negative rows
298    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
299      double value = scalar*x[iColumn];
300      if (value) {
301        int iRowM = indices_[j];
302        int iRowP = indices_[j+1];
303        if (iRowM>=0)
304          y[iRowM] -= value;
305        if (iRowP>=0)
306          y[iRowP] += value;
307      }
308    }
309  }
310}
311void 
312ClpNetworkMatrix::transposeTimes(double scalar,
313                                const double * x, double * y) const
314{
315  int iColumn;
316  CoinBigIndex j=0;
317  if (trueNetwork_) {
318    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
319      double value = y[iColumn];
320      int iRowM = indices_[j];
321      int iRowP = indices_[j+1];
322      value -= scalar*x[iRowM];
323      value += scalar*x[iRowP];
324      y[iColumn] = value;
325    }
326  } else {
327    // skip negative rows
328    for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
329      double value = y[iColumn];
330      int iRowM = indices_[j];
331      int iRowP = indices_[j+1];
332      if (iRowM>=0)
333        value -= scalar*x[iRowM];
334      if (iRowP>=0)
335        value += scalar*x[iRowP];
336      y[iColumn] = value;
337    }
338  }
339}
340void 
341ClpNetworkMatrix::times(double scalar,
342                       const double * x, double * y,
343                       const double * rowScale, 
344                       const double * columnScale) const
345{
346  // we know it is not scaled
347  times(scalar, x, y);
348}
349void 
350ClpNetworkMatrix::transposeTimes( double scalar,
351                                 const double * x, double * y,
352                                 const double * rowScale, 
353                                 const double * columnScale, double * spare) const
354{
355  // we know it is not scaled
356  transposeTimes(scalar, x, y);
357}
358/* Return <code>x * A + y</code> in <code>z</code>.
359        Squashes small elements and knows about ClpSimplex */
360void 
361ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
362                              const CoinIndexedVector * rowArray,
363                              CoinIndexedVector * y,
364                              CoinIndexedVector * columnArray) const
365{
366  // we know it is not scaled
367  columnArray->clear();
368  double * pi = rowArray->denseVector();
369  int numberNonZero=0;
370  int * index = columnArray->getIndices();
371  double * array = columnArray->denseVector();
372  int numberInRowArray = rowArray->getNumElements();
373  // maybe I need one in OsiSimplex
374  double zeroTolerance = model->factorization()->zeroTolerance();
375  int numberRows = model->numberRows();
376#ifndef NO_RTTI
377  ClpPlusMinusOneMatrix* rowCopy =
378    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
379#else
380  ClpPlusMinusOneMatrix* rowCopy =
381    static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
382#endif
383  bool packed = rowArray->packedMode();
384  double factor = 0.3;
385  // We may not want to do by row if there may be cache problems
386  int numberColumns = model->numberColumns();
387  // It would be nice to find L2 cache size - for moment 512K
388  // Be slightly optimistic
389  if (numberColumns*sizeof(double)>1000000) {
390    if (numberRows*10<numberColumns)
391      factor=0.1;
392    else if (numberRows*4<numberColumns)
393      factor=0.15;
394    else if (numberRows*2<numberColumns)
395      factor=0.2;
396    //if (model->numberIterations()%50==0)
397    //printf("%d nonzero\n",numberInRowArray);
398  }
399  if (numberInRowArray>factor*numberRows||!rowCopy) {
400    // do by column
401    int iColumn;
402    assert (!y->getNumElements());
403    CoinBigIndex j=0;
404    if (packed) {
405      // need to expand pi into y
406      assert(y->capacity()>=numberRows);
407      double * piOld = pi;
408      pi = y->denseVector();
409      const int * whichRow = rowArray->getIndices();
410      int i;
411      // modify pi so can collapse to one loop
412      for (i=0;i<numberInRowArray;i++) {
413        int iRow = whichRow[i];
414        pi[iRow]=scalar*piOld[i];
415      }
416      if (trueNetwork_) {
417        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
418          double value = 0.0;
419          int iRowM = indices_[j];
420          int iRowP = indices_[j+1];
421          value -= pi[iRowM];
422          value += pi[iRowP];
423          if (fabs(value)>zeroTolerance) {
424            array[numberNonZero]=value;
425            index[numberNonZero++]=iColumn;
426          }
427        }
428      } else {
429        // skip negative rows
430        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
431          double value = 0.0;
432          int iRowM = indices_[j];
433          int iRowP = indices_[j+1];
434          if (iRowM>=0)
435            value -= pi[iRowM];
436          if (iRowP>=0)
437            value += pi[iRowP];
438          if (fabs(value)>zeroTolerance) {
439            array[numberNonZero]=value;
440            index[numberNonZero++]=iColumn;
441          }
442        }
443      }
444      for (i=0;i<numberInRowArray;i++) {
445        int iRow = whichRow[i];
446        pi[iRow]=0.0;
447      }
448    } else {
449      if (trueNetwork_) {
450        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
451          double value = 0.0;
452          int iRowM = indices_[j];
453          int iRowP = indices_[j+1];
454          value -= scalar*pi[iRowM];
455          value += scalar*pi[iRowP];
456          if (fabs(value)>zeroTolerance) {
457            index[numberNonZero++]=iColumn;
458            array[iColumn]=value;
459          }
460        }
461      } else {
462        // skip negative rows
463        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
464          double value = 0.0;
465          int iRowM = indices_[j];
466          int iRowP = indices_[j+1];
467          if (iRowM>=0)
468            value -= scalar*pi[iRowM];
469          if (iRowP>=0)
470            value += scalar*pi[iRowP];
471          if (fabs(value)>zeroTolerance) {
472            index[numberNonZero++]=iColumn;
473            array[iColumn]=value;
474          }
475        }
476      }
477    }
478    columnArray->setNumElements(numberNonZero);
479  } else {
480    // do by row
481    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
482  }
483}
484/* Return <code>x *A in <code>z</code> but
485   just for indices in y. */
486void 
487ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * model,
488                              const CoinIndexedVector * rowArray,
489                              const CoinIndexedVector * y,
490                              CoinIndexedVector * columnArray) const
491{
492  columnArray->clear();
493  double * pi = rowArray->denseVector();
494  double * array = columnArray->denseVector();
495  int jColumn;
496  int numberToDo = y->getNumElements();
497  const int * which = y->getIndices();
498  assert (!rowArray->packedMode());
499  columnArray->setPacked();
500  if (trueNetwork_) {
501    for (jColumn=0;jColumn<numberToDo;jColumn++) {
502      int iColumn = which[jColumn];
503      double value = 0.0;
504      CoinBigIndex j=iColumn<<1;
505      int iRowM = indices_[j];
506      int iRowP = indices_[j+1];
507      value -= pi[iRowM];
508      value += pi[iRowP];
509      array[jColumn]=value;
510    }
511  } else {
512    // skip negative rows
513    for (jColumn=0;jColumn<numberToDo;jColumn++) {
514      int iColumn = which[jColumn];
515      double value = 0.0;
516      CoinBigIndex j=iColumn<<1;
517      int iRowM = indices_[j];
518      int iRowP = indices_[j+1];
519      if (iRowM>=0)
520        value -= pi[iRowM];
521      if (iRowP>=0)
522        value += pi[iRowP];
523      array[jColumn]=value;
524    }
525  }
526}
527/// returns number of elements in column part of basis,
528CoinBigIndex
529ClpNetworkMatrix::countBasis(ClpSimplex * model,
530                                 const int * whichColumn, 
531                                 int numberBasic,
532                             int & numberColumnBasic)
533{
534  int i;
535  CoinBigIndex numberElements=0;
536  if (trueNetwork_) {
537    numberElements = 2*numberColumnBasic;
538  } else {
539    for (i=0;i<numberColumnBasic;i++) {
540      int iColumn = whichColumn[i];
541      CoinBigIndex j=iColumn<<1;
542      int iRowM = indices_[j];
543      int iRowP = indices_[j+1];
544      if (iRowM>=0)
545        numberElements ++;
546      if (iRowP>=0)
547        numberElements ++;
548    }
549  }
550  return numberElements;
551}
552void
553ClpNetworkMatrix::fillBasis(ClpSimplex * model,
554                         const int * whichColumn, 
555                         int & numberColumnBasic,
556                         int * indexRowU, int * start,
557                         int * rowCount, int * columnCount,
558                         double * elementU)
559{
560  int i;
561  CoinBigIndex numberElements=start[0];
562  if (trueNetwork_) {
563    for (i=0;i<numberColumnBasic;i++) {
564      int iColumn = whichColumn[i];
565      CoinBigIndex j=iColumn<<1;
566      int iRowM = indices_[j];
567      int iRowP = indices_[j+1];
568      indexRowU[numberElements]=iRowM;
569      rowCount[iRowM++];
570      elementU[numberElements]=-1.0;
571      indexRowU[numberElements+1]=iRowP;
572      rowCount[iRowP++];
573      elementU[numberElements+1]=1.0;
574      numberElements+=2;
575      start[i+1]=numberElements;
576      columnCount[i]=2;
577    }
578  } else {
579    for (i=0;i<numberColumnBasic;i++) {
580      int iColumn = whichColumn[i];
581      CoinBigIndex j=iColumn<<1;
582      int iRowM = indices_[j];
583      int iRowP = indices_[j+1];
584      if (iRowM>=0) {
585        indexRowU[numberElements]=iRowM;
586        rowCount[iRowM++];
587        elementU[numberElements++]=-1.0;
588      }
589      if (iRowP>=0) {
590        indexRowU[numberElements]=iRowP;
591        rowCount[iRowP++];
592        elementU[numberElements++]=1.0;
593      }
594      start[i+1]=numberElements;
595      columnCount[i]=numberElements-start[i];
596    }
597  }
598}
599/* Unpacks a column into an CoinIndexedvector
600 */
601void 
602ClpNetworkMatrix::unpack(const ClpSimplex * model,CoinIndexedVector * rowArray,
603                   int iColumn) const 
604{
605  CoinBigIndex j=iColumn<<1;
606  int iRowM = indices_[j];
607  int iRowP = indices_[j+1];
608  if (iRowM>=0) 
609    rowArray->add(iRowM,-1.0);
610  if (iRowP>=0) 
611    rowArray->add(iRowP,1.0);
612}
613/* Unpacks a column into an CoinIndexedvector
614** in packed foramt
615Note that model is NOT const.  Bounds and objective could
616be modified if doing column generation (just for this variable) */
617void 
618ClpNetworkMatrix::unpackPacked(ClpSimplex * model,
619                            CoinIndexedVector * rowArray,
620                            int iColumn) const
621{
622  int * index = rowArray->getIndices();
623  double * array = rowArray->denseVector();
624  int number = 0;
625  CoinBigIndex j=iColumn<<1;
626  int iRowM = indices_[j];
627  int iRowP = indices_[j+1];
628  if (iRowM>=0) {
629    array[number]=-1.0;
630    index[number++]=iRowM;
631  }
632  if (iRowP>=0) {
633    array[number]=1.0;
634    index[number++]=iRowP;
635  }
636  rowArray->setNumElements(number);
637  rowArray->setPackedMode(true);
638}
639/* Adds multiple of a column into an CoinIndexedvector
640      You can use quickAdd to add to vector */
641void 
642ClpNetworkMatrix::add(const ClpSimplex * model,CoinIndexedVector * rowArray,
643                   int iColumn, double multiplier) const 
644{
645  CoinBigIndex j=iColumn<<1;
646  int iRowM = indices_[j];
647  int iRowP = indices_[j+1];
648  if (iRowM>=0) 
649    rowArray->quickAdd(iRowM,-multiplier);
650  if (iRowP>=0) 
651    rowArray->quickAdd(iRowP,multiplier);
652}
653/* Adds multiple of a column into an array */
654void 
655ClpNetworkMatrix::add(const ClpSimplex * model,double * array,
656                    int iColumn, double multiplier) const
657{
658  CoinBigIndex j=iColumn<<1;
659  int iRowM = indices_[j];
660  int iRowP = indices_[j+1];
661  if (iRowM>=0) 
662    array[iRowM] -= multiplier;
663  if (iRowP>=0) 
664    array[iRowP] += multiplier;
665}
666
667// Return a complete CoinPackedMatrix
668CoinPackedMatrix * 
669ClpNetworkMatrix::getPackedMatrix() const 
670{
671  return new CoinPackedMatrix(true,numberRows_,numberColumns_,
672                              2*numberColumns_,
673                              getElements(),indices_,
674                              getVectorStarts(),getVectorLengths());
675
676}
677/* A vector containing the elements in the packed matrix. Note that there
678   might be gaps in this list, entries that do not belong to any
679   major-dimension vector. To get the actual elements one should look at
680   this vector together with vectorStarts and vectorLengths. */
681const double * 
682ClpNetworkMatrix::getElements() const 
683{
684  assert (trueNetwork_); // fix later
685  if (!elements_) {
686    elements_ = new double [2*numberColumns_];
687    int i;
688    for (i=0;i<2*numberColumns_;i+=2) {
689      elements_[i]=-1.0;
690      elements_[i+1]=1.0;
691    }
692  }
693  return elements_;
694}
695
696const CoinBigIndex * 
697ClpNetworkMatrix::getVectorStarts() const 
698{
699  assert (trueNetwork_); // fix later
700  if (!starts_) {
701    starts_ = new CoinBigIndex [numberColumns_+1];
702    int i;
703    for (i=0;i<numberColumns_+1;i++) {
704      starts_[i]=2*i;
705    }
706  }
707  return starts_;
708}
709/* The lengths of the major-dimension vectors. */
710const int * 
711ClpNetworkMatrix::getVectorLengths() const
712{
713  assert (trueNetwork_); // fix later
714  if (!lengths_) {
715    lengths_ = new int [numberColumns_];
716    int i;
717    for (i=0;i<numberColumns_;i++) {
718      lengths_[i]=2;
719    }
720  }
721  return lengths_;
722}
723/* Delete the columns whose indices are listed in <code>indDel</code>. */
724void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel) 
725{
726  std::cerr<<"deleteCols not implemented in ClpNetworkMatrix"<<std::endl;
727  abort();
728}
729/* Delete the rows whose indices are listed in <code>indDel</code>. */
730void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel) 
731{
732  std::cerr<<"deleteRows not implemented in ClpNetworkMatrix"<<std::endl;
733  abort();
734}
735/* Given positive integer weights for each row fills in sum of weights
736   for each column (and slack).
737   Returns weights vector
738*/
739CoinBigIndex * 
740ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
741{
742  int numberRows = model->numberRows();
743  int numberColumns =model->numberColumns();
744  int number = numberRows+numberColumns;
745  CoinBigIndex * weights = new CoinBigIndex[number];
746  int i;
747  for (i=0;i<numberColumns;i++) {
748    CoinBigIndex j=i<<1;
749    CoinBigIndex count=0;
750    int iRowM = indices_[j];
751    int iRowP = indices_[j+1];
752    if (iRowM>=0) {
753      count += inputWeights[iRowM];
754    }
755    if (iRowP>=0) {
756      count += inputWeights[iRowP];
757    }
758    weights[i]=count;
759  }
760  for (i=0;i<numberRows;i++) {
761    weights[i+numberColumns]=inputWeights[i];
762  }
763  return weights;
764}
765/* Returns largest and smallest elements of both signs.
766   Largest refers to largest absolute value.
767*/
768void 
769ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
770                       double & smallestPositive, double & largestPositive)
771{
772  smallestNegative=-1.0;
773  largestNegative=-1.0;
774  smallestPositive=1.0;
775  largestPositive=1.0;
776}
777// Says whether it can do partial pricing
778bool 
779ClpNetworkMatrix::canDoPartialPricing() const
780{
781  return true; 
782}
783// Partial pricing
784void 
785ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
786                              int & bestSequence, int & numberWanted)
787{
788  numberWanted=currentWanted_;
789  int j;
790  int start = (int) (startFraction*numberColumns_);
791  int end = CoinMin((int) (endFraction*numberColumns_+1),numberColumns_);
792  double tolerance=model->currentDualTolerance();
793  double * reducedCost = model->djRegion();
794  const double * duals = model->dualRowSolution();
795  const double * cost = model->costRegion();
796  double bestDj;
797  if (bestSequence>=0)
798    bestDj = fabs(reducedCost[bestSequence]);
799  else
800    bestDj=tolerance;
801  int sequenceOut = model->sequenceOut();
802  int saveSequence = bestSequence;
803  if (!trueNetwork_) {
804    // Not true network
805    int iSequence;
806    for (iSequence=start;iSequence<end;iSequence++) {
807      if (iSequence!=sequenceOut) {
808        double value;
809        int iRowM,iRowP;
810        ClpSimplex::Status status = model->getStatus(iSequence);
811       
812        switch(status) {
813         
814        case ClpSimplex::basic:
815        case ClpSimplex::isFixed:
816          break;
817        case ClpSimplex::isFree:
818        case ClpSimplex::superBasic:
819          value=cost[iSequence];
820          j = iSequence<<1;
821          // skip negative rows
822          iRowM = indices_[j];
823          iRowP = indices_[j+1];
824          if (iRowM>=0)
825            value += duals[iRowM];
826          if (iRowP>=0)
827            value -= duals[iRowP];
828          value = fabs(value);
829          if (value>FREE_ACCEPT*tolerance) {
830            numberWanted--;
831            // we are going to bias towards free (but only if reasonable)
832            value *= FREE_BIAS;
833            if (value>bestDj) {
834              // check flagged variable and correct dj
835              if (!model->flagged(iSequence)) {
836                bestDj=value;
837                bestSequence = iSequence;
838              } else {
839                // just to make sure we don't exit before got something
840                numberWanted++;
841              }
842            }
843          }
844          break;
845        case ClpSimplex::atUpperBound:
846          value=cost[iSequence];
847          j = iSequence<<1;
848          // skip negative rows
849          iRowM = indices_[j];
850          iRowP = indices_[j+1];
851          if (iRowM>=0)
852            value += duals[iRowM];
853          if (iRowP>=0)
854            value -= duals[iRowP];
855          if (value>tolerance) {
856            numberWanted--;
857            if (value>bestDj) {
858              // check flagged variable and correct dj
859              if (!model->flagged(iSequence)) {
860                bestDj=value;
861                bestSequence = iSequence;
862              } else {
863                // just to make sure we don't exit before got something
864                numberWanted++;
865              }
866            }
867          }
868          break;
869        case ClpSimplex::atLowerBound:
870          value=cost[iSequence];
871          j = iSequence<<1;
872          // skip negative rows
873          iRowM = indices_[j];
874          iRowP = indices_[j+1];
875          if (iRowM>=0)
876            value += duals[iRowM];
877          if (iRowP>=0)
878            value -= duals[iRowP];
879          value = -value;
880          if (value>tolerance) {
881            numberWanted--;
882            if (value>bestDj) {
883              // check flagged variable and correct dj
884              if (!model->flagged(iSequence)) {
885                bestDj=value;
886                bestSequence = iSequence;
887              } else {
888                // just to make sure we don't exit before got something
889                numberWanted++;
890              }
891            }
892          }
893          break;
894        }
895      }
896      if (!numberWanted)
897        break;
898    }
899    if (bestSequence!=saveSequence) {
900      // recompute dj
901      double value=cost[bestSequence];
902      j = bestSequence<<1;
903      // skip negative rows
904      int iRowM = indices_[j];
905      int iRowP = indices_[j+1];
906      if (iRowM>=0)
907        value += duals[iRowM];
908      if (iRowP>=0)
909        value -= duals[iRowP];
910      reducedCost[bestSequence] = value;
911      savedBestSequence_ = bestSequence;
912      savedBestDj_ = reducedCost[savedBestSequence_];
913    }
914  } else {
915    // true network
916    int iSequence;
917    for (iSequence=start;iSequence<end;iSequence++) {
918      if (iSequence!=sequenceOut) {
919        double value;
920        int iRowM,iRowP;
921        ClpSimplex::Status status = model->getStatus(iSequence);
922       
923        switch(status) {
924         
925        case ClpSimplex::basic:
926        case ClpSimplex::isFixed:
927          break;
928        case ClpSimplex::isFree:
929        case ClpSimplex::superBasic:
930          value=cost[iSequence];
931          j = iSequence<<1;
932          iRowM = indices_[j];
933          iRowP = indices_[j+1];
934          value += duals[iRowM];
935          value -= duals[iRowP];
936          value = fabs(value);
937          if (value>FREE_ACCEPT*tolerance) {
938            numberWanted--;
939            // we are going to bias towards free (but only if reasonable)
940            value *= FREE_BIAS;
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::atUpperBound:
954          value=cost[iSequence];
955          j = iSequence<<1;
956          iRowM = indices_[j];
957          iRowP = indices_[j+1];
958          value += duals[iRowM];
959          value -= duals[iRowP];
960          if (value>tolerance) {
961            numberWanted--;
962            if (value>bestDj) {
963              // check flagged variable and correct dj
964              if (!model->flagged(iSequence)) {
965                bestDj=value;
966                bestSequence = iSequence;
967              } else {
968                // just to make sure we don't exit before got something
969                numberWanted++;
970              }
971            }
972          }
973          break;
974        case ClpSimplex::atLowerBound:
975          value=cost[iSequence];
976          j = iSequence<<1;
977          iRowM = indices_[j];
978          iRowP = indices_[j+1];
979          value += duals[iRowM];
980          value -= duals[iRowP];
981          value = -value;
982          if (value>tolerance) {
983            numberWanted--;
984            if (value>bestDj) {
985              // check flagged variable and correct dj
986              if (!model->flagged(iSequence)) {
987                bestDj=value;
988                bestSequence = iSequence;
989              } else {
990                // just to make sure we don't exit before got something
991                numberWanted++;
992              }
993            }
994          }
995          break;
996        }
997      }
998      if (!numberWanted)
999        break;
1000    }
1001    if (bestSequence!=saveSequence) {
1002      // recompute dj
1003      double value=cost[bestSequence];
1004      j = bestSequence<<1;
1005      int iRowM = indices_[j];
1006      int iRowP = indices_[j+1];
1007      value += duals[iRowM];
1008      value -= duals[iRowP];
1009      reducedCost[bestSequence] = value;
1010      savedBestSequence_ = bestSequence;
1011      savedBestDj_ = reducedCost[savedBestSequence_];
1012    }
1013  }
1014  currentWanted_=numberWanted;
1015}
1016// Allow any parts of a created CoinMatrix to be deleted
1017void 
1018ClpNetworkMatrix::releasePackedMatrix() const 
1019{
1020  delete [] elements_;
1021  delete [] lengths_;
1022  elements_=NULL;
1023  lengths_=NULL;
1024}
Note: See TracBrowser for help on using the repository browser.