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

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

get rid of compiler warnings

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 33.2 KB
Line 
1/* $Id: ClpNetworkMatrix.cpp 1402 2009-07-25 08:39:55Z 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*/, 
350                                  double * /*spare*/) const
351{
352  // we know it is not scaled
353  transposeTimes(scalar, x, y);
354}
355/* Return <code>x * A + y</code> in <code>z</code>.
356        Squashes small elements and knows about ClpSimplex */
357void 
358ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
359                              const CoinIndexedVector * rowArray,
360                              CoinIndexedVector * y,
361                              CoinIndexedVector * columnArray) const
362{
363  // we know it is not scaled
364  columnArray->clear();
365  double * pi = rowArray->denseVector();
366  int numberNonZero=0;
367  int * index = columnArray->getIndices();
368  double * array = columnArray->denseVector();
369  int numberInRowArray = rowArray->getNumElements();
370  // maybe I need one in OsiSimplex
371  double zeroTolerance = model->zeroTolerance();
372  int numberRows = model->numberRows();
373#ifndef NO_RTTI
374  ClpPlusMinusOneMatrix* rowCopy =
375    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
376#else
377  ClpPlusMinusOneMatrix* rowCopy =
378    static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
379#endif
380  bool packed = rowArray->packedMode();
381  double factor = 0.3;
382  // We may not want to do by row if there may be cache problems
383  int numberColumns = model->numberColumns();
384  // It would be nice to find L2 cache size - for moment 512K
385  // Be slightly optimistic
386  if (numberColumns*sizeof(double)>1000000) {
387    if (numberRows*10<numberColumns)
388      factor=0.1;
389    else if (numberRows*4<numberColumns)
390      factor=0.15;
391    else if (numberRows*2<numberColumns)
392      factor=0.2;
393    //if (model->numberIterations()%50==0)
394    //printf("%d nonzero\n",numberInRowArray);
395  }
396  if (numberInRowArray>factor*numberRows||!rowCopy) {
397    // do by column
398    int iColumn;
399    assert (!y->getNumElements());
400    CoinBigIndex j=0;
401    if (packed) {
402      // need to expand pi into y
403      assert(y->capacity()>=numberRows);
404      double * piOld = pi;
405      pi = y->denseVector();
406      const int * whichRow = rowArray->getIndices();
407      int i;
408      // modify pi so can collapse to one loop
409      for (i=0;i<numberInRowArray;i++) {
410        int iRow = whichRow[i];
411        pi[iRow]=scalar*piOld[i];
412      }
413      if (trueNetwork_) {
414        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
415          double value = 0.0;
416          int iRowM = indices_[j];
417          int iRowP = indices_[j+1];
418          value -= pi[iRowM];
419          value += pi[iRowP];
420          if (fabs(value)>zeroTolerance) {
421            array[numberNonZero]=value;
422            index[numberNonZero++]=iColumn;
423          }
424        }
425      } else {
426        // skip negative rows
427        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
428          double value = 0.0;
429          int iRowM = indices_[j];
430          int iRowP = indices_[j+1];
431          if (iRowM>=0)
432            value -= pi[iRowM];
433          if (iRowP>=0)
434            value += pi[iRowP];
435          if (fabs(value)>zeroTolerance) {
436            array[numberNonZero]=value;
437            index[numberNonZero++]=iColumn;
438          }
439        }
440      }
441      for (i=0;i<numberInRowArray;i++) {
442        int iRow = whichRow[i];
443        pi[iRow]=0.0;
444      }
445    } else {
446      if (trueNetwork_) {
447        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
448          double value = 0.0;
449          int iRowM = indices_[j];
450          int iRowP = indices_[j+1];
451          value -= scalar*pi[iRowM];
452          value += scalar*pi[iRowP];
453          if (fabs(value)>zeroTolerance) {
454            index[numberNonZero++]=iColumn;
455            array[iColumn]=value;
456          }
457        }
458      } else {
459        // skip negative rows
460        for (iColumn=0;iColumn<numberColumns_;iColumn++,j+=2) {
461          double value = 0.0;
462          int iRowM = indices_[j];
463          int iRowP = indices_[j+1];
464          if (iRowM>=0)
465            value -= scalar*pi[iRowM];
466          if (iRowP>=0)
467            value += scalar*pi[iRowP];
468          if (fabs(value)>zeroTolerance) {
469            index[numberNonZero++]=iColumn;
470            array[iColumn]=value;
471          }
472        }
473      }
474    }
475    columnArray->setNumElements(numberNonZero);
476  } else {
477    // do by row
478    rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
479  }
480}
481/* Return <code>x *A in <code>z</code> but
482   just for indices in y. */
483void 
484ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/,
485                              const CoinIndexedVector * rowArray,
486                              const CoinIndexedVector * y,
487                              CoinIndexedVector * columnArray) const
488{
489  columnArray->clear();
490  double * pi = rowArray->denseVector();
491  double * array = columnArray->denseVector();
492  int jColumn;
493  int numberToDo = y->getNumElements();
494  const int * which = y->getIndices();
495  assert (!rowArray->packedMode());
496  columnArray->setPacked();
497  if (trueNetwork_) {
498    for (jColumn=0;jColumn<numberToDo;jColumn++) {
499      int iColumn = which[jColumn];
500      double value = 0.0;
501      CoinBigIndex j=iColumn<<1;
502      int iRowM = indices_[j];
503      int iRowP = indices_[j+1];
504      value -= pi[iRowM];
505      value += pi[iRowP];
506      array[jColumn]=value;
507    }
508  } else {
509    // skip negative rows
510    for (jColumn=0;jColumn<numberToDo;jColumn++) {
511      int iColumn = which[jColumn];
512      double value = 0.0;
513      CoinBigIndex j=iColumn<<1;
514      int iRowM = indices_[j];
515      int iRowP = indices_[j+1];
516      if (iRowM>=0)
517        value -= pi[iRowM];
518      if (iRowP>=0)
519        value += pi[iRowP];
520      array[jColumn]=value;
521    }
522  }
523}
524/// returns number of elements in column part of basis,
525CoinBigIndex
526ClpNetworkMatrix::countBasis( const int * whichColumn, 
527                             int & numberColumnBasic)
528{
529  int i;
530  CoinBigIndex numberElements=0;
531  if (trueNetwork_) {
532    numberElements = 2*numberColumnBasic;
533  } else {
534    for (i=0;i<numberColumnBasic;i++) {
535      int iColumn = whichColumn[i];
536      CoinBigIndex j=iColumn<<1;
537      int iRowM = indices_[j];
538      int iRowP = indices_[j+1];
539      if (iRowM>=0)
540        numberElements ++;
541      if (iRowP>=0)
542        numberElements ++;
543    }
544  }
545  return numberElements;
546}
547void
548ClpNetworkMatrix::fillBasis(ClpSimplex * /*model*/,
549                         const int * whichColumn, 
550                         int & numberColumnBasic,
551                         int * indexRowU, int * start,
552                         int * rowCount, int * columnCount,
553                         CoinFactorizationDouble * elementU)
554{
555  int i;
556  CoinBigIndex numberElements=start[0];
557  if (trueNetwork_) {
558    for (i=0;i<numberColumnBasic;i++) {
559      int iColumn = whichColumn[i];
560      CoinBigIndex j=iColumn<<1;
561      int iRowM = indices_[j];
562      int iRowP = indices_[j+1];
563      indexRowU[numberElements]=iRowM;
564      rowCount[iRowM]++;
565      elementU[numberElements]=-1.0;
566      indexRowU[numberElements+1]=iRowP;
567      rowCount[iRowP]++;
568      elementU[numberElements+1]=1.0;
569      numberElements+=2;
570      start[i+1]=numberElements;
571      columnCount[i]=2;
572    }
573  } else {
574    for (i=0;i<numberColumnBasic;i++) {
575      int iColumn = whichColumn[i];
576      CoinBigIndex j=iColumn<<1;
577      int iRowM = indices_[j];
578      int iRowP = indices_[j+1];
579      if (iRowM>=0) {
580        indexRowU[numberElements]=iRowM;
581        rowCount[iRowM]++;
582        elementU[numberElements++]=-1.0;
583      }
584      if (iRowP>=0) {
585        indexRowU[numberElements]=iRowP;
586        rowCount[iRowP]++;
587        elementU[numberElements++]=1.0;
588      }
589      start[i+1]=numberElements;
590      columnCount[i]=numberElements-start[i];
591    }
592  }
593}
594/* Unpacks a column into an CoinIndexedvector
595 */
596void 
597ClpNetworkMatrix::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 
613ClpNetworkMatrix::unpackPacked(ClpSimplex * /*model*/,
614                            CoinIndexedVector * rowArray,
615                            int iColumn) const
616{
617  int * index = rowArray->getIndices();
618  double * array = rowArray->denseVector();
619  int number = 0;
620  CoinBigIndex j=iColumn<<1;
621  int iRowM = indices_[j];
622  int iRowP = indices_[j+1];
623  if (iRowM>=0) {
624    array[number]=-1.0;
625    index[number++]=iRowM;
626  }
627  if (iRowP>=0) {
628    array[number]=1.0;
629    index[number++]=iRowP;
630  }
631  rowArray->setNumElements(number);
632  rowArray->setPackedMode(true);
633}
634/* Adds multiple of a column into an CoinIndexedvector
635      You can use quickAdd to add to vector */
636void 
637ClpNetworkMatrix::add(const ClpSimplex * /*model*/,CoinIndexedVector * rowArray,
638                   int iColumn, double multiplier) const 
639{
640  CoinBigIndex j=iColumn<<1;
641  int iRowM = indices_[j];
642  int iRowP = indices_[j+1];
643  if (iRowM>=0) 
644    rowArray->quickAdd(iRowM,-multiplier);
645  if (iRowP>=0) 
646    rowArray->quickAdd(iRowP,multiplier);
647}
648/* Adds multiple of a column into an array */
649void 
650ClpNetworkMatrix::add(const ClpSimplex * /*model*/,double * array,
651                    int iColumn, double multiplier) const
652{
653  CoinBigIndex j=iColumn<<1;
654  int iRowM = indices_[j];
655  int iRowP = indices_[j+1];
656  if (iRowM>=0) 
657    array[iRowM] -= multiplier;
658  if (iRowP>=0) 
659    array[iRowP] += multiplier;
660}
661
662// Return a complete CoinPackedMatrix
663CoinPackedMatrix * 
664ClpNetworkMatrix::getPackedMatrix() const 
665{
666  if (!matrix_) {
667    assert (trueNetwork_); // fix later
668    int numberElements = 2*numberColumns_;
669    double * elements = new double [numberElements];
670    CoinBigIndex i;
671    for (i=0;i<2*numberColumns_;i+=2) {
672      elements[i]=-1.0;
673      elements[i+1]=1.0;
674    }
675    CoinBigIndex * starts = new CoinBigIndex [numberColumns_+1];
676    for (i=0;i<numberColumns_+1;i++) {
677      starts[i]=2*i;
678    }
679    // use assignMatrix to save space
680    delete [] lengths_;
681    lengths_ = NULL;
682    matrix_ = new CoinPackedMatrix();
683    int * indices = CoinCopyOfArray(indices_,2*numberColumns_);
684    matrix_->assignMatrix(true,numberRows_,numberColumns_,
685                          getNumElements(),
686                          elements,indices,
687                          starts,lengths_);
688    assert(!elements);
689    assert(!starts);
690    assert (!indices);
691    assert (!lengths_);
692  }
693  return matrix_;
694}
695/* A vector containing the elements in the packed matrix. Note that there
696   might be gaps in this list, entries that do not belong to any
697   major-dimension vector. To get the actual elements one should look at
698   this vector together with vectorStarts and vectorLengths. */
699const double * 
700ClpNetworkMatrix::getElements() const 
701{
702  if (!matrix_) 
703    getPackedMatrix();
704  return matrix_->getElements();
705}
706
707const CoinBigIndex * 
708ClpNetworkMatrix::getVectorStarts() const 
709{
710  if (!matrix_) 
711    getPackedMatrix();
712  return matrix_->getVectorStarts();
713}
714/* The lengths of the major-dimension vectors. */
715const int * 
716ClpNetworkMatrix::getVectorLengths() const
717{
718  assert (trueNetwork_); // fix later
719  if (!lengths_) {
720    lengths_ = new int [numberColumns_];
721    int i;
722    for (i=0;i<numberColumns_;i++) {
723      lengths_[i]=2;
724    }
725  }
726  return lengths_;
727}
728/* Delete the columns whose indices are listed in <code>indDel</code>. */
729void ClpNetworkMatrix::deleteCols(const int numDel, const int * indDel) 
730{
731  assert (trueNetwork_);
732  int iColumn;
733  int numberBad=0;
734  // Use array to make sure we can have duplicates
735  char * which = new char[numberColumns_];
736  memset(which,0,numberColumns_);
737  int nDuplicate=0;
738  for (iColumn=0;iColumn<numDel;iColumn++) {
739    int jColumn = indDel[iColumn];
740    if (jColumn<0||jColumn>=numberColumns_) {
741      numberBad++;
742    } else {
743      if (which[jColumn])
744        nDuplicate++;
745      else
746        which[jColumn]=1;
747    }
748  }
749  if (numberBad)
750    throw CoinError("Indices out of range", "deleteCols", "ClpNetworkMatrix");
751  int newNumber = numberColumns_-numDel+nDuplicate;
752  // Get rid of temporary arrays
753  delete [] lengths_;
754  lengths_=NULL;
755  delete matrix_;
756  matrix_= NULL;
757  int newSize=2*newNumber;
758  int * newIndices = new int [newSize];
759  newSize=0;
760  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
761    if (!which[iColumn]) {
762      CoinBigIndex start;
763      CoinBigIndex i;
764      start = 2*iColumn;
765      for (i=start;i<start+2;i++) 
766        newIndices[newSize++]=indices_[i];
767    }
768  }
769  delete [] which;
770  delete [] indices_;
771  indices_= newIndices;
772  numberColumns_ = newNumber;
773}
774/* Delete the rows whose indices are listed in <code>indDel</code>. */
775void ClpNetworkMatrix::deleteRows(const int numDel, const int * indDel) 
776{
777  int iRow;
778  int numberBad=0;
779  // Use array to make sure we can have duplicates
780  int * which = new int [numberRows_];
781  memset(which,0,numberRows_*sizeof(int));
782  for (iRow=0;iRow<numDel;iRow++) {
783    int jRow = indDel[iRow];
784    if (jRow<0||jRow>=numberRows_) {
785      numberBad++;
786    } else {
787      which[jRow]=1;
788    }
789  }
790  if (numberBad)
791    throw CoinError("Indices out of range", "deleteRows", "ClpNetworkMatrix");
792  // Only valid of all columns have 0 entries
793  int iColumn;
794  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
795    CoinBigIndex start;
796    CoinBigIndex i;
797    start = 2*iColumn;
798    for (i=start;i<start+2;i++) {
799      int iRow = indices_[i];
800      if (which[iRow])
801        numberBad++;
802    }
803  }
804  if (numberBad)
805    throw CoinError("Row has entries", "deleteRows", "ClpNetworkMatrix");
806  int newNumber=0;
807  for (iRow=0;iRow<numberRows_;iRow++) {
808    if (!which[iRow])
809      which[iRow]=newNumber++;
810    else
811      which[iRow]=-1;
812  }
813  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
814    CoinBigIndex start;
815    CoinBigIndex i;
816    start = 2*iColumn;
817    for (i=start;i<start+2;i++) {
818      int iRow = indices_[i];
819      indices_[i]=which[iRow];
820    }
821  }
822  delete [] which;
823  numberRows_ = newNumber;
824}
825/* Given positive integer weights for each row fills in sum of weights
826   for each column (and slack).
827   Returns weights vector
828*/
829CoinBigIndex * 
830ClpNetworkMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
831{
832  int numberRows = model->numberRows();
833  int numberColumns =model->numberColumns();
834  int number = numberRows+numberColumns;
835  CoinBigIndex * weights = new CoinBigIndex[number];
836  int i;
837  for (i=0;i<numberColumns;i++) {
838    CoinBigIndex j=i<<1;
839    CoinBigIndex count=0;
840    int iRowM = indices_[j];
841    int iRowP = indices_[j+1];
842    if (iRowM>=0) {
843      count += inputWeights[iRowM];
844    }
845    if (iRowP>=0) {
846      count += inputWeights[iRowP];
847    }
848    weights[i]=count;
849  }
850  for (i=0;i<numberRows;i++) {
851    weights[i+numberColumns]=inputWeights[i];
852  }
853  return weights;
854}
855/* Returns largest and smallest elements of both signs.
856   Largest refers to largest absolute value.
857*/
858void 
859ClpNetworkMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
860                       double & smallestPositive, double & largestPositive)
861{
862  smallestNegative=-1.0;
863  largestNegative=-1.0;
864  smallestPositive=1.0;
865  largestPositive=1.0;
866}
867// Says whether it can do partial pricing
868bool 
869ClpNetworkMatrix::canDoPartialPricing() const
870{
871  return true; 
872}
873// Partial pricing
874void 
875ClpNetworkMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
876                              int & bestSequence, int & numberWanted)
877{
878  numberWanted=currentWanted_;
879  int j;
880  int start = static_cast<int> (startFraction*numberColumns_);
881  int end = CoinMin(static_cast<int> (endFraction*numberColumns_+1),numberColumns_);
882  double tolerance=model->currentDualTolerance();
883  double * reducedCost = model->djRegion();
884  const double * duals = model->dualRowSolution();
885  const double * cost = model->costRegion();
886  double bestDj;
887  if (bestSequence>=0)
888    bestDj = fabs(reducedCost[bestSequence]);
889  else
890    bestDj=tolerance;
891  int sequenceOut = model->sequenceOut();
892  int saveSequence = bestSequence;
893  if (!trueNetwork_) {
894    // Not true network
895    int iSequence;
896    for (iSequence=start;iSequence<end;iSequence++) {
897      if (iSequence!=sequenceOut) {
898        double value;
899        int iRowM,iRowP;
900        ClpSimplex::Status status = model->getStatus(iSequence);
901
902        switch(status) {
903         
904        case ClpSimplex::basic:
905        case ClpSimplex::isFixed:
906          break;
907        case ClpSimplex::isFree:
908        case ClpSimplex::superBasic:
909          value=cost[iSequence];
910          j = iSequence<<1;
911          // skip negative rows
912          iRowM = indices_[j];
913          iRowP = indices_[j+1];
914          if (iRowM>=0)
915            value += duals[iRowM];
916          if (iRowP>=0)
917            value -= duals[iRowP];
918          value = fabs(value);
919          if (value>FREE_ACCEPT*tolerance) {
920            numberWanted--;
921            // we are going to bias towards free (but only if reasonable)
922            value *= FREE_BIAS;
923            if (value>bestDj) {
924              // check flagged variable and correct dj
925              if (!model->flagged(iSequence)) {
926                bestDj=value;
927                bestSequence = iSequence;
928              } else {
929                // just to make sure we don't exit before got something
930                numberWanted++;
931              }
932            }
933          }
934          break;
935        case ClpSimplex::atUpperBound:
936          value=cost[iSequence];
937          j = iSequence<<1;
938          // skip negative rows
939          iRowM = indices_[j];
940          iRowP = indices_[j+1];
941          if (iRowM>=0)
942            value += duals[iRowM];
943          if (iRowP>=0)
944            value -= duals[iRowP];
945          if (value>tolerance) {
946            numberWanted--;
947            if (value>bestDj) {
948              // check flagged variable and correct dj
949              if (!model->flagged(iSequence)) {
950                bestDj=value;
951                bestSequence = iSequence;
952              } else {
953                // just to make sure we don't exit before got something
954                numberWanted++;
955              }
956            }
957          }
958          break;
959        case ClpSimplex::atLowerBound:
960          value=cost[iSequence];
961          j = iSequence<<1;
962          // skip negative rows
963          iRowM = indices_[j];
964          iRowP = indices_[j+1];
965          if (iRowM>=0)
966            value += duals[iRowM];
967          if (iRowP>=0)
968            value -= duals[iRowP];
969          value = -value;
970          if (value>tolerance) {
971            numberWanted--;
972            if (value>bestDj) {
973              // check flagged variable and correct dj
974              if (!model->flagged(iSequence)) {
975                bestDj=value;
976                bestSequence = iSequence;
977              } else {
978                // just to make sure we don't exit before got something
979                numberWanted++;
980              }
981            }
982          }
983          break;
984        }
985      }
986      if (!numberWanted)
987        break;
988    }
989    if (bestSequence!=saveSequence) {
990      // recompute dj
991      double value=cost[bestSequence];
992      j = bestSequence<<1;
993      // skip negative rows
994      int iRowM = indices_[j];
995      int iRowP = indices_[j+1];
996      if (iRowM>=0)
997        value += duals[iRowM];
998      if (iRowP>=0)
999        value -= duals[iRowP];
1000      reducedCost[bestSequence] = value;
1001      savedBestSequence_ = bestSequence;
1002      savedBestDj_ = reducedCost[savedBestSequence_];
1003    }
1004  } else {
1005    // true network
1006    int iSequence;
1007    for (iSequence=start;iSequence<end;iSequence++) {
1008      if (iSequence!=sequenceOut) {
1009        double value;
1010        int iRowM,iRowP;
1011        ClpSimplex::Status status = model->getStatus(iSequence);
1012
1013        switch(status) {
1014         
1015        case ClpSimplex::basic:
1016        case ClpSimplex::isFixed:
1017          break;
1018        case ClpSimplex::isFree:
1019        case ClpSimplex::superBasic:
1020          value=cost[iSequence];
1021          j = iSequence<<1;
1022          iRowM = indices_[j];
1023          iRowP = indices_[j+1];
1024          value += duals[iRowM];
1025          value -= duals[iRowP];
1026          value = fabs(value);
1027          if (value>FREE_ACCEPT*tolerance) {
1028            numberWanted--;
1029            // we are going to bias towards free (but only if reasonable)
1030            value *= FREE_BIAS;
1031            if (value>bestDj) {
1032              // check flagged variable and correct dj
1033              if (!model->flagged(iSequence)) {
1034                bestDj=value;
1035                bestSequence = iSequence;
1036              } else {
1037                // just to make sure we don't exit before got something
1038                numberWanted++;
1039              }
1040            }
1041          }
1042          break;
1043        case ClpSimplex::atUpperBound:
1044          value=cost[iSequence];
1045          j = iSequence<<1;
1046          iRowM = indices_[j];
1047          iRowP = indices_[j+1];
1048          value += duals[iRowM];
1049          value -= duals[iRowP];
1050          if (value>tolerance) {
1051            numberWanted--;
1052            if (value>bestDj) {
1053              // check flagged variable and correct dj
1054              if (!model->flagged(iSequence)) {
1055                bestDj=value;
1056                bestSequence = iSequence;
1057              } else {
1058                // just to make sure we don't exit before got something
1059                numberWanted++;
1060              }
1061            }
1062          }
1063          break;
1064        case ClpSimplex::atLowerBound:
1065          value=cost[iSequence];
1066          j = iSequence<<1;
1067          iRowM = indices_[j];
1068          iRowP = indices_[j+1];
1069          value += duals[iRowM];
1070          value -= duals[iRowP];
1071          value = -value;
1072          if (value>tolerance) {
1073            numberWanted--;
1074            if (value>bestDj) {
1075              // check flagged variable and correct dj
1076              if (!model->flagged(iSequence)) {
1077                bestDj=value;
1078                bestSequence = iSequence;
1079              } else {
1080                // just to make sure we don't exit before got something
1081                numberWanted++;
1082              }
1083            }
1084          }
1085          break;
1086        }
1087      }
1088      if (!numberWanted)
1089        break;
1090    }
1091    if (bestSequence!=saveSequence) {
1092      // recompute dj
1093      double value=cost[bestSequence];
1094      j = bestSequence<<1;
1095      int iRowM = indices_[j];
1096      int iRowP = indices_[j+1];
1097      value += duals[iRowM];
1098      value -= duals[iRowP];
1099      reducedCost[bestSequence] = value;
1100      savedBestSequence_ = bestSequence;
1101      savedBestDj_ = reducedCost[savedBestSequence_];
1102    }
1103  }
1104  currentWanted_=numberWanted;
1105}
1106// Allow any parts of a created CoinMatrix to be deleted
1107void 
1108ClpNetworkMatrix::releasePackedMatrix() const 
1109{
1110  delete matrix_;
1111  delete [] lengths_;
1112  matrix_=NULL;
1113  lengths_=NULL;
1114}
1115// Append Columns
1116void 
1117ClpNetworkMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
1118{
1119  int iColumn;
1120  int numberBad=0;
1121  for (iColumn=0;iColumn<number;iColumn++) {
1122    int n=columns[iColumn]->getNumElements();
1123    const double * element = columns[iColumn]->getElements();
1124    if (n!=2)
1125      numberBad++;
1126    if (fabs(element[0])!=1.0||fabs(element[1])!=1.0)
1127        numberBad++;
1128    else if (element[0]*element[1]!=-1.0)
1129        numberBad++;
1130  }
1131  if (numberBad)
1132    throw CoinError("Not network", "appendCols", "ClpNetworkMatrix");
1133  // Get rid of temporary arrays
1134  delete [] lengths_;
1135  lengths_=NULL;
1136  delete matrix_;
1137  matrix_= NULL;
1138  CoinBigIndex size = 2*number;
1139  int * temp2 = new int [numberColumns_*2+size];
1140  CoinMemcpyN(indices_,numberColumns_*2,temp2);
1141  delete [] indices_;
1142  indices_= temp2;
1143  // now add
1144  size=2*numberColumns_;
1145  for (iColumn=0;iColumn<number;iColumn++) {
1146    const int * row = columns[iColumn]->getIndices();
1147    const double * element = columns[iColumn]->getElements();
1148    if (element[0]==-1.0) {
1149      indices_[size++]=row[0];
1150      indices_[size++]=row[1];
1151    } else {
1152      indices_[size++]=row[1];
1153      indices_[size++]=row[0];
1154    }
1155  }
1156 
1157  numberColumns_ += number;
1158}
1159// Append Rows
1160void 
1161ClpNetworkMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
1162{
1163  // must be zero arrays
1164  int numberBad=0;
1165  int iRow;
1166  for (iRow=0;iRow<number;iRow++) {
1167    numberBad += rows[iRow]->getNumElements();
1168  }
1169  if (numberBad)
1170    throw CoinError("Not NULL rows", "appendRows", "ClpNetworkMatrix");
1171  numberRows_ += number;
1172}
1173#ifndef SLIM_CLP
1174/* Append a set of rows/columns to the end of the matrix. Returns number of errors
1175   i.e. if any of the new rows/columns contain an index that's larger than the
1176   number of columns-1/rows-1 (if numberOther>0) or duplicates
1177   If 0 then rows, 1 if columns */
1178int 
1179ClpNetworkMatrix::appendMatrix(int number, int type,
1180                                    const CoinBigIndex * starts, const int * index,
1181                               const double * element, int /*numberOther*/)
1182{
1183  int numberErrors=0;
1184  // make into CoinPackedVector
1185  CoinPackedVectorBase ** vectors=
1186    new CoinPackedVectorBase * [number];
1187  int iVector;
1188  for (iVector=0;iVector<number;iVector++) {
1189    int iStart = starts[iVector];
1190    vectors[iVector] = 
1191      new CoinPackedVector(starts[iVector+1]-iStart,
1192                           index+iStart,element+iStart);
1193  }
1194  if (type==0) {
1195    // rows
1196    appendRows(number,vectors);
1197  } else {
1198    // columns
1199    appendCols(number,vectors);
1200  }
1201  for (iVector=0;iVector<number;iVector++) 
1202    delete vectors[iVector];
1203  delete [] vectors;
1204  return numberErrors;
1205}
1206#endif
1207/* Subset clone (without gaps).  Duplicates are allowed
1208   and order is as given */
1209ClpMatrixBase * 
1210ClpNetworkMatrix::subsetClone (int numberRows, const int * whichRows,
1211                              int numberColumns, 
1212                              const int * whichColumns) const 
1213{
1214  return new ClpNetworkMatrix(*this, numberRows, whichRows,
1215                                   numberColumns, whichColumns);
1216}
1217/* Subset constructor (without gaps).  Duplicates are allowed
1218   and order is as given */
1219ClpNetworkMatrix::ClpNetworkMatrix (
1220                       const ClpNetworkMatrix & rhs,
1221                       int numberRows, const int * whichRow,
1222                       int numberColumns, const int * whichColumn)
1223: ClpMatrixBase(rhs)
1224{ 
1225  setType(11);
1226  matrix_ = NULL;
1227  lengths_=NULL;
1228  indices_=new int[2*numberColumns];;
1229  numberRows_=numberRows;
1230  numberColumns_=numberColumns;
1231  trueNetwork_=true;
1232  int iColumn;
1233  int numberBad=0;
1234  int * which = new int [rhs.numberRows_];
1235  int iRow;
1236  for (iRow=0;iRow<rhs.numberRows_;iRow++) 
1237    which[iRow]=-1;
1238  int n=0;
1239  for (iRow=0;iRow<numberRows;iRow++) {
1240    int jRow=whichRow[iRow];
1241    assert (jRow>=0&&jRow<rhs.numberRows_);
1242    which[jRow]=n++;
1243  }
1244  for (iColumn=0;iColumn<numberColumns;iColumn++) {
1245    CoinBigIndex start;
1246    CoinBigIndex i;
1247    start = 2*iColumn;
1248    CoinBigIndex offset = 2*whichColumn[iColumn]-start;
1249    for (i=start;i<start+2;i++) {
1250      int iRow = rhs.indices_[i+offset];
1251      iRow = which[iRow];
1252      if (iRow<0)
1253        numberBad++;
1254      else
1255        indices_[i]=iRow;
1256    }
1257  }
1258  if (numberBad)
1259    throw CoinError("Invalid rows", "subsetConstructor", "ClpNetworkMatrix");
1260}
Note: See TracBrowser for help on using the repository browser.