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

Last change on this file since 2030 was 1665, checked in by lou, 9 years ago

Add EPL license notice in src.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 45.0 KB
Line 
1/* $Id: ClpNetworkMatrix.cpp 1665 2011-01-04 17:55:54Z forrest $ */
2// Copyright (C) 2003, International Business Machines
3// Corporation and others.  All Rights Reserved.
4// This code is licensed under the terms of the Eclipse Public License (EPL).
5
6
7#include <cstdio>
8
9#include "CoinPragma.hpp"
10#include "CoinIndexedVector.hpp"
11#include "CoinHelperFunctions.hpp"
12#include "CoinPackedVector.hpp"
13
14#include "ClpSimplex.hpp"
15#include "ClpFactorization.hpp"
16// at end to get min/max!
17#include "ClpNetworkMatrix.hpp"
18#include "ClpPlusMinusOneMatrix.hpp"
19#include "ClpMessage.hpp"
20#include <iostream>
21#include <cassert>
22
23//#############################################################################
24// Constructors / Destructor / Assignment
25//#############################################################################
26
27//-------------------------------------------------------------------
28// Default Constructor
29//-------------------------------------------------------------------
30ClpNetworkMatrix::ClpNetworkMatrix ()
31     : ClpMatrixBase()
32{
33     setType(11);
34     matrix_ = NULL;
35     lengths_ = NULL;
36     indices_ = NULL;
37     numberRows_ = 0;
38     numberColumns_ = 0;
39     trueNetwork_ = false;
40}
41
42/* Constructor from two arrays */
43ClpNetworkMatrix::ClpNetworkMatrix(int numberColumns, const int * head,
44                                   const int * tail)
45     : ClpMatrixBase()
46{
47     setType(11);
48     matrix_ = NULL;
49     lengths_ = NULL;
50     indices_ = new int[2*numberColumns];;
51     numberRows_ = -1;
52     numberColumns_ = numberColumns;
53     trueNetwork_ = true;
54     int iColumn;
55     CoinBigIndex j = 0;
56     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
57          int iRow = head[iColumn];
58          numberRows_ = CoinMax(numberRows_, iRow);
59          indices_[j] = iRow;
60          iRow = tail[iColumn];
61          numberRows_ = CoinMax(numberRows_, iRow);
62          indices_[j+1] = iRow;
63     }
64     numberRows_++;
65}
66//-------------------------------------------------------------------
67// Copy constructor
68//-------------------------------------------------------------------
69ClpNetworkMatrix::ClpNetworkMatrix (const ClpNetworkMatrix & rhs)
70     : ClpMatrixBase(rhs)
71{
72     matrix_ = NULL;
73     lengths_ = NULL;
74     indices_ = NULL;
75     numberRows_ = rhs.numberRows_;
76     numberColumns_ = rhs.numberColumns_;
77     trueNetwork_ = rhs.trueNetwork_;
78     if (numberColumns_) {
79          indices_ = new int [ 2*numberColumns_];
80          CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
81     }
82     int numberRows = getNumRows();
83     if (rhs.rhsOffset_ && numberRows) {
84          rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
85     } else {
86          rhsOffset_ = NULL;
87     }
88}
89
90ClpNetworkMatrix::ClpNetworkMatrix (const CoinPackedMatrix & rhs)
91     : ClpMatrixBase()
92{
93     setType(11);
94     matrix_ = NULL;
95     lengths_ = NULL;
96     indices_ = NULL;
97     int iColumn;
98     assert (rhs.isColOrdered());
99     // get matrix data pointers
100     const int * row = rhs.getIndices();
101     const CoinBigIndex * columnStart = rhs.getVectorStarts();
102     const int * columnLength = rhs.getVectorLengths();
103     const double * elementByColumn = rhs.getElements();
104     numberColumns_ = rhs.getNumCols();
105     int goodNetwork = 1;
106     numberRows_ = -1;
107     indices_ = new int[2*numberColumns_];
108     CoinBigIndex j = 0;
109     for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
110          CoinBigIndex k = columnStart[iColumn];
111          int iRow;
112          switch (columnLength[iColumn]) {
113          case 0:
114               goodNetwork = -1; // not classic network
115               indices_[j] = -1;
116               indices_[j+1] = -1;
117               break;
118
119          case 1:
120               goodNetwork = -1; // not classic network
121               if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
122                    indices_[j] = -1;
123                    iRow = row[k];
124                    numberRows_ = CoinMax(numberRows_, iRow);
125                    indices_[j+1] = iRow;
126               } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
127                    indices_[j+1] = -1;
128                    iRow = row[k];
129                    numberRows_ = CoinMax(numberRows_, iRow);
130                    indices_[j] = iRow;
131               } else {
132                    goodNetwork = 0; // not a network
133               }
134               break;
135
136          case 2:
137               if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
138                    if (fabs(elementByColumn[k+1] + 1.0) < 1.0e-10) {
139                         iRow = row[k];
140                         numberRows_ = CoinMax(numberRows_, iRow);
141                         indices_[j+1] = iRow;
142                         iRow = row[k+1];
143                         numberRows_ = CoinMax(numberRows_, iRow);
144                         indices_[j] = iRow;
145                    } else {
146                         goodNetwork = 0; // not a network
147                    }
148               } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
149                    if (fabs(elementByColumn[k+1] - 1.0) < 1.0e-10) {
150                         iRow = row[k];
151                         numberRows_ = CoinMax(numberRows_, iRow);
152                         indices_[j] = iRow;
153                         iRow = row[k+1];
154                         numberRows_ = CoinMax(numberRows_, iRow);
155                         indices_[j+1] = iRow;
156                    } else {
157                         goodNetwork = 0; // not a network
158                    }
159               } else {
160                    goodNetwork = 0; // not a network
161               }
162               break;
163
164          default:
165               goodNetwork = 0; // not a network
166               break;
167          }
168          if (!goodNetwork)
169               break;
170     }
171     if (!goodNetwork) {
172          delete [] indices_;
173          // put in message
174          printf("Not a network - can test if indices_ null\n");
175          indices_ = NULL;
176          numberRows_ = 0;
177          numberColumns_ = 0;
178     } else {
179          numberRows_ ++; //  correct
180          trueNetwork_ = goodNetwork > 0;
181     }
182}
183
184//-------------------------------------------------------------------
185// Destructor
186//-------------------------------------------------------------------
187ClpNetworkMatrix::~ClpNetworkMatrix ()
188{
189     delete matrix_;
190     delete [] lengths_;
191     delete [] indices_;
192}
193
194//----------------------------------------------------------------
195// Assignment operator
196//-------------------------------------------------------------------
197ClpNetworkMatrix &
198ClpNetworkMatrix::operator=(const ClpNetworkMatrix& rhs)
199{
200     if (this != &rhs) {
201          ClpMatrixBase::operator=(rhs);
202          delete matrix_;
203          delete [] lengths_;
204          delete [] indices_;
205          matrix_ = NULL;
206          lengths_ = NULL;
207          indices_ = NULL;
208          numberRows_ = rhs.numberRows_;
209          numberColumns_ = rhs.numberColumns_;
210          trueNetwork_ = rhs.trueNetwork_;
211          if (numberColumns_) {
212               indices_ = new int [ 2*numberColumns_];
213               CoinMemcpyN(rhs.indices_, 2 * numberColumns_, indices_);
214          }
215     }
216     return *this;
217}
218//-------------------------------------------------------------------
219// Clone
220//-------------------------------------------------------------------
221ClpMatrixBase * ClpNetworkMatrix::clone() const
222{
223     return new ClpNetworkMatrix(*this);
224}
225
226/* Returns a new matrix in reverse order without gaps */
227ClpMatrixBase *
228ClpNetworkMatrix::reverseOrderedCopy() const
229{
230     // count number in each row
231     CoinBigIndex * tempP = new CoinBigIndex [numberRows_];
232     CoinBigIndex * tempN = new CoinBigIndex [numberRows_];
233     memset(tempP, 0, numberRows_ * sizeof(CoinBigIndex));
234     memset(tempN, 0, numberRows_ * sizeof(CoinBigIndex));
235     CoinBigIndex j = 0;
236     int i;
237     for (i = 0; i < numberColumns_; i++, j += 2) {
238          int iRow = indices_[j];
239          tempN[iRow]++;
240          iRow = indices_[j+1];
241          tempP[iRow]++;
242     }
243     int * newIndices = new int [2*numberColumns_];
244     CoinBigIndex * newP = new CoinBigIndex [numberRows_+1];
245     CoinBigIndex * newN = new CoinBigIndex[numberRows_];
246     int iRow;
247     j = 0;
248     // do starts
249     for (iRow = 0; iRow < numberRows_; iRow++) {
250          newP[iRow] = j;
251          j += tempP[iRow];
252          tempP[iRow] = newP[iRow];
253          newN[iRow] = j;
254          j += tempN[iRow];
255          tempN[iRow] = newN[iRow];
256     }
257     newP[numberRows_] = j;
258     j = 0;
259     for (i = 0; i < numberColumns_; i++, j += 2) {
260          int iRow = indices_[j];
261          CoinBigIndex put = tempN[iRow];
262          newIndices[put++] = i;
263          tempN[iRow] = put;
264          iRow = indices_[j+1];
265          put = tempP[iRow];
266          newIndices[put++] = i;
267          tempP[iRow] = put;
268     }
269     delete [] tempP;
270     delete [] tempN;
271     ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
272     newCopy->passInCopy(numberRows_, numberColumns_,
273                         false,  newIndices, newP, newN);
274     return newCopy;
275}
276//unscaled versions
277void
278ClpNetworkMatrix::times(double scalar,
279                        const double * x, double * y) const
280{
281     int iColumn;
282     CoinBigIndex j = 0;
283     if (trueNetwork_) {
284          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
285               double value = scalar * x[iColumn];
286               if (value) {
287                    int iRowM = indices_[j];
288                    int iRowP = indices_[j+1];
289                    y[iRowM] -= value;
290                    y[iRowP] += value;
291               }
292          }
293     } else {
294          // skip negative rows
295          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
296               double value = scalar * x[iColumn];
297               if (value) {
298                    int iRowM = indices_[j];
299                    int iRowP = indices_[j+1];
300                    if (iRowM >= 0)
301                         y[iRowM] -= value;
302                    if (iRowP >= 0)
303                         y[iRowP] += value;
304               }
305          }
306     }
307}
308void
309ClpNetworkMatrix::transposeTimes(double scalar,
310                                 const double * x, double * y) const
311{
312     int iColumn;
313     CoinBigIndex j = 0;
314     if (trueNetwork_) {
315          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
316               double value = y[iColumn];
317               int iRowM = indices_[j];
318               int iRowP = indices_[j+1];
319               value -= scalar * x[iRowM];
320               value += scalar * x[iRowP];
321               y[iColumn] = value;
322          }
323     } else {
324          // skip negative rows
325          for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
326               double value = y[iColumn];
327               int iRowM = indices_[j];
328               int iRowP = indices_[j+1];
329               if (iRowM >= 0)
330                    value -= scalar * x[iRowM];
331               if (iRowP >= 0)
332                    value += scalar * x[iRowP];
333               y[iColumn] = value;
334          }
335     }
336}
337void
338ClpNetworkMatrix::times(double scalar,
339                        const double * x, double * y,
340                        const double * /*rowScale*/,
341                        const double * /*columnScale*/) const
342{
343     // we know it is not scaled
344     times(scalar, x, y);
345}
346void
347ClpNetworkMatrix::transposeTimes( double scalar,
348                                  const double * x, double * y,
349                                  const double * /*rowScale*/,
350                                  const double * /*columnScale*/,
351                                  double * /*spare*/) const
352{
353     // we know it is not scaled
354     transposeTimes(scalar, x, y);
355}
356/* Return <code>x * A + y</code> in <code>z</code>.
357        Squashes small elements and knows about ClpSimplex */
358void
359ClpNetworkMatrix::transposeTimes(const ClpSimplex * model, double scalar,
360                                 const CoinIndexedVector * rowArray,
361                                 CoinIndexedVector * y,
362                                 CoinIndexedVector * columnArray) const
363{
364     // we know it is not scaled
365     columnArray->clear();
366     double * pi = rowArray->denseVector();
367     int numberNonZero = 0;
368     int * index = columnArray->getIndices();
369     double * array = columnArray->denseVector();
370     int numberInRowArray = rowArray->getNumElements();
371     // maybe I need one in OsiSimplex
372     double zeroTolerance = model->zeroTolerance();
373     int numberRows = model->numberRows();
374#ifndef NO_RTTI
375     ClpPlusMinusOneMatrix* rowCopy =
376          dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
377#else
378     ClpPlusMinusOneMatrix* rowCopy =
379          static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
380#endif
381     bool packed = rowArray->packedMode();
382     double factor = 0.3;
383     // We may not want to do by row if there may be cache problems
384     int numberColumns = model->numberColumns();
385     // It would be nice to find L2 cache size - for moment 512K
386     // Be slightly optimistic
387     if (numberColumns * sizeof(double) > 1000000) {
388          if (numberRows * 10 < numberColumns)
389               factor = 0.1;
390          else if (numberRows * 4 < numberColumns)
391               factor = 0.15;
392          else if (numberRows * 2 < numberColumns)
393               factor = 0.2;
394          //if (model->numberIterations()%50==0)
395          //printf("%d nonzero\n",numberInRowArray);
396     }
397     if (numberInRowArray > factor * numberRows || !rowCopy) {
398          // do by column
399          int iColumn;
400          assert (!y->getNumElements());
401          CoinBigIndex j = 0;
402          if (packed) {
403               // need to expand pi into y
404               assert(y->capacity() >= numberRows);
405               double * piOld = pi;
406               pi = y->denseVector();
407               const int * whichRow = rowArray->getIndices();
408               int i;
409               // modify pi so can collapse to one loop
410               for (i = 0; i < numberInRowArray; i++) {
411                    int iRow = whichRow[i];
412                    pi[iRow] = scalar * piOld[i];
413               }
414               if (trueNetwork_) {
415                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
416                         double value = 0.0;
417                         int iRowM = indices_[j];
418                         int iRowP = indices_[j+1];
419                         value -= pi[iRowM];
420                         value += pi[iRowP];
421                         if (fabs(value) > zeroTolerance) {
422                              array[numberNonZero] = value;
423                              index[numberNonZero++] = iColumn;
424                         }
425                    }
426               } else {
427                    // skip negative rows
428                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
429                         double value = 0.0;
430                         int iRowM = indices_[j];
431                         int iRowP = indices_[j+1];
432                         if (iRowM >= 0)
433                              value -= pi[iRowM];
434                         if (iRowP >= 0)
435                              value += pi[iRowP];
436                         if (fabs(value) > zeroTolerance) {
437                              array[numberNonZero] = value;
438                              index[numberNonZero++] = iColumn;
439                         }
440                    }
441               }
442               for (i = 0; i < numberInRowArray; i++) {
443                    int iRow = whichRow[i];
444                    pi[iRow] = 0.0;
445               }
446          } else {
447               if (trueNetwork_) {
448                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
449                         double value = 0.0;
450                         int iRowM = indices_[j];
451                         int iRowP = indices_[j+1];
452                         value -= scalar * pi[iRowM];
453                         value += scalar * pi[iRowP];
454                         if (fabs(value) > zeroTolerance) {
455                              index[numberNonZero++] = iColumn;
456                              array[iColumn] = value;
457                         }
458                    }
459               } else {
460                    // skip negative rows
461                    for (iColumn = 0; iColumn < numberColumns_; iColumn++, j += 2) {
462                         double value = 0.0;
463                         int iRowM = indices_[j];
464                         int iRowP = indices_[j+1];
465                         if (iRowM >= 0)
466                              value -= scalar * pi[iRowM];
467                         if (iRowP >= 0)
468                              value += scalar * pi[iRowP];
469                         if (fabs(value) > zeroTolerance) {
470                              index[numberNonZero++] = iColumn;
471                              array[iColumn] = value;
472                         }
473                    }
474               }
475          }
476          columnArray->setNumElements(numberNonZero);
477     } else {
478          // do by row
479          rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
480     }
481}
482/* Return <code>x *A in <code>z</code> but
483   just for indices in y. */
484void
485ClpNetworkMatrix::subsetTransposeTimes(const ClpSimplex * /*model*/,
486                                       const CoinIndexedVector * rowArray,
487                                       const CoinIndexedVector * y,
488                                       CoinIndexedVector * columnArray) const
489{
490     columnArray->clear();
491     double * pi = rowArray->denseVector();
492     double * array = columnArray->denseVector();
493     int jColumn;
494     int numberToDo = y->getNumElements();
495     const int * which = y->getIndices();
496     assert (!rowArray->packedMode());
497     columnArray->setPacked();
498     if (trueNetwork_) {
499          for (jColumn = 0; jColumn < numberToDo; jColumn++) {
500               int iColumn = which[jColumn];
501               double value = 0.0;
502               CoinBigIndex j = iColumn << 1;
503               int iRowM = indices_[j];
504               int iRowP = indices_[j+1];
505               value -= pi[iRowM];
506               value += pi[iRowP];
507               array[jColumn] = value;
508          }
509     } else {
510          // skip negative rows
511          for (jColumn = 0; jColumn < numberToDo; jColumn++) {
512               int iColumn = which[jColumn];
513               double value = 0.0;
514               CoinBigIndex j = iColumn << 1;
515               int iRowM = indices_[j];
516               int iRowP = indices_[j+1];
517               if (iRowM >= 0)
518                    value -= pi[iRowM];
519               if (iRowP >= 0)
520                    value += pi[iRowP];
521               array[jColumn] = value;
522          }
523     }
524}
525/// returns number of elements in column part of basis,
526CoinBigIndex
527ClpNetworkMatrix::countBasis( const int * whichColumn,
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.