source: trunk/Clp/src/ClpPlusMinusOneMatrix.cpp @ 2271

Last change on this file since 2271 was 2271, checked in by forrest, 2 years ago

change some ints to CoinBigIndex?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Id
File size: 121.9 KB
Line 
1/* $Id: ClpPlusMinusOneMatrix.cpp 2271 2017-08-14 08:50:39Z 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 "CoinPackedVector.hpp"
12#include "CoinHelperFunctions.hpp"
13
14#include "ClpSimplex.hpp"
15#include "ClpFactorization.hpp"
16// at end to get min/max!
17#include "ClpPlusMinusOneMatrix.hpp"
18#include "ClpMessage.hpp"
19#ifdef CLP_PLUS_ONE_MATRIX
20static int oneitcount[13]={0,0,0,0,0,0,0,0,0,0,0,0,0};
21static void oneit(int i)
22{
23  if (!oneitcount[i])
24    printf("Plus ones for call %d\n",i);
25  oneitcount[i]++;
26  oneitcount[12]++;
27  if ((oneitcount[12]%1000)==0) {
28    printf("Plus counts");
29    for (int j=0;j<12;j++)
30      printf(" %d",oneitcount[j]);
31    printf("\n");
32  }   
33}
34#endif
35//#############################################################################
36// Constructors / Destructor / Assignment
37//#############################################################################
38
39//-------------------------------------------------------------------
40// Default Constructor
41//-------------------------------------------------------------------
42ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix ()
43     : ClpMatrixBase()
44{
45     setType(12);
46     matrix_ = NULL;
47     startPositive_ = NULL;
48     startNegative_ = NULL;
49     lengths_ = NULL;
50     indices_ = NULL;
51     numberRows_ = 0;
52     numberColumns_ = 0;
53#ifdef CLP_PLUS_ONE_MATRIX
54     otherFlags_ = 0;
55#endif
56     columnOrdered_ = true;
57}
58
59//-------------------------------------------------------------------
60// Copy constructor
61//-------------------------------------------------------------------
62ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const ClpPlusMinusOneMatrix & rhs)
63     : ClpMatrixBase(rhs)
64{
65     matrix_ = NULL;
66     startPositive_ = NULL;
67     startNegative_ = NULL;
68     lengths_ = NULL;
69     indices_ = NULL;
70     numberRows_ = rhs.numberRows_;
71     numberColumns_ = rhs.numberColumns_;
72#ifdef CLP_PLUS_ONE_MATRIX
73     otherFlags_ = rhs.otherFlags_;
74#endif
75     columnOrdered_ = rhs.columnOrdered_;
76     if (numberColumns_) {
77          CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
78          indices_ = new int [ numberElements];
79          CoinMemcpyN(rhs.indices_, numberElements, indices_);
80          startPositive_ = new CoinBigIndex [ numberColumns_+1];
81          CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_);
82          startNegative_ = new CoinBigIndex [ numberColumns_];
83          CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_);
84     }
85     int numberRows = getNumRows();
86     if (rhs.rhsOffset_ && numberRows) {
87          rhsOffset_ = ClpCopyOfArray(rhs.rhsOffset_, numberRows);
88     } else {
89          rhsOffset_ = NULL;
90     }
91}
92// Constructor from arrays
93ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix(int numberRows, int numberColumns,
94          bool columnOrdered, const int * indices,
95          const CoinBigIndex * startPositive,
96          const CoinBigIndex * startNegative)
97     : ClpMatrixBase()
98{
99     setType(12);
100     matrix_ = NULL;
101     lengths_ = NULL;
102     numberRows_ = numberRows;
103     numberColumns_ = numberColumns;
104     columnOrdered_ = columnOrdered;
105#ifdef CLP_PLUS_ONE_MATRIX
106     otherFlags_ = 0;
107#endif
108     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
109     CoinBigIndex numberElements = startPositive[numberMajor];
110     startPositive_ = ClpCopyOfArray(startPositive, numberMajor + 1);
111     startNegative_ = ClpCopyOfArray(startNegative, numberMajor);
112     indices_ = ClpCopyOfArray(indices, numberElements);
113     // Check valid
114     checkValid(false);
115}
116
117ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (const CoinPackedMatrix & rhs)
118     : ClpMatrixBase()
119{
120     setType(12);
121     matrix_ = NULL;
122     startPositive_ = NULL;
123     startNegative_ = NULL;
124     lengths_ = NULL;
125     indices_ = NULL;
126     int iColumn;
127     assert (rhs.isColOrdered());
128     // get matrix data pointers
129     const int * row = rhs.getIndices();
130     const CoinBigIndex * columnStart = rhs.getVectorStarts();
131     const int * columnLength = rhs.getVectorLengths();
132     const double * elementByColumn = rhs.getElements();
133     numberColumns_ = rhs.getNumCols();
134     numberRows_ = -1;
135     indices_ = new int[rhs.getNumElements()];
136     startPositive_ = new CoinBigIndex [numberColumns_+1];
137     startNegative_ = new CoinBigIndex [numberColumns_];
138     int * temp = new int [rhs.getNumRows()];
139     CoinBigIndex j = 0;
140     CoinBigIndex numberGoodP = 0;
141     CoinBigIndex numberGoodM = 0;
142     CoinBigIndex numberBad = 0;
143     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
144          CoinBigIndex k;
145          int iNeg = 0;
146          startPositive_[iColumn] = j;
147          for (k = columnStart[iColumn]; k < columnStart[iColumn] + columnLength[iColumn];
148                    k++) {
149               int iRow;
150               if (fabs(elementByColumn[k] - 1.0) < 1.0e-10) {
151                    iRow = row[k];
152                    numberRows_ = CoinMax(numberRows_, iRow);
153                    indices_[j++] = iRow;
154                    numberGoodP++;
155               } else if (fabs(elementByColumn[k] + 1.0) < 1.0e-10) {
156                    iRow = row[k];
157                    numberRows_ = CoinMax(numberRows_, iRow);
158                    temp[iNeg++] = iRow;
159                    numberGoodM++;
160               } else {
161                    numberBad++;
162               }
163          }
164          // move negative
165          startNegative_[iColumn] = j;
166          for (k = 0; k < iNeg; k++) {
167               indices_[j++] = temp[k];
168          }
169     }
170     startPositive_[numberColumns_] = j;
171     delete [] temp;
172     if (numberBad) {
173          delete [] indices_;
174          indices_ = NULL;
175          numberRows_ = 0;
176          numberColumns_ = 0;
177          delete [] startPositive_;
178          delete [] startNegative_;
179          // Put in statistics
180          startPositive_ = new CoinBigIndex [3];
181          startPositive_[0] = numberGoodP;
182          startPositive_[1] = numberGoodM;
183          startPositive_[2] = numberBad;
184          startNegative_ = NULL;
185     } else {
186          numberRows_ ++; //  correct
187          // but number should be same as rhs
188          assert (numberRows_ <= rhs.getNumRows());
189          numberRows_ = rhs.getNumRows();
190          columnOrdered_ = true;
191     }
192     // Check valid
193     if (!numberBad)
194          checkValid(false);
195}
196
197//-------------------------------------------------------------------
198// Destructor
199//-------------------------------------------------------------------
200ClpPlusMinusOneMatrix::~ClpPlusMinusOneMatrix ()
201{
202     delete matrix_;
203     delete [] startPositive_;
204     delete [] startNegative_;
205     delete [] lengths_;
206     delete [] indices_;
207}
208
209//----------------------------------------------------------------
210// Assignment operator
211//-------------------------------------------------------------------
212ClpPlusMinusOneMatrix &
213ClpPlusMinusOneMatrix::operator=(const ClpPlusMinusOneMatrix& rhs)
214{
215     if (this != &rhs) {
216          ClpMatrixBase::operator=(rhs);
217          delete matrix_;
218          delete [] startPositive_;
219          delete [] startNegative_;
220          delete [] lengths_;
221          delete [] indices_;
222          matrix_ = NULL;
223          startPositive_ = NULL;
224          lengths_ = NULL;
225          indices_ = NULL;
226          numberRows_ = rhs.numberRows_;
227          numberColumns_ = rhs.numberColumns_;
228          columnOrdered_ = rhs.columnOrdered_;
229#ifdef CLP_PLUS_ONE_MATRIX
230          otherFlags_ = rhs.otherFlags_;
231#endif
232          if (numberColumns_) {
233               CoinBigIndex numberElements = rhs.startPositive_[numberColumns_];
234               indices_ = new int [ numberElements];
235               CoinMemcpyN(rhs.indices_, numberElements, indices_);
236               startPositive_ = new CoinBigIndex [ numberColumns_+1];
237               CoinMemcpyN(rhs.startPositive_, (numberColumns_ + 1), startPositive_);
238               startNegative_ = new CoinBigIndex [ numberColumns_];
239               CoinMemcpyN(rhs.startNegative_, numberColumns_, startNegative_);
240          }
241     }
242     return *this;
243}
244//-------------------------------------------------------------------
245// Clone
246//-------------------------------------------------------------------
247ClpMatrixBase * ClpPlusMinusOneMatrix::clone() const
248{
249     return new ClpPlusMinusOneMatrix(*this);
250}
251/* Subset clone (without gaps).  Duplicates are allowed
252   and order is as given */
253ClpMatrixBase *
254ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
255                                    int numberColumns,
256                                    const int * whichColumns) const
257{
258     return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
259                                      numberColumns, whichColumns);
260}
261/* Subset constructor (without gaps).  Duplicates are allowed
262   and order is as given */
263ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
264     const ClpPlusMinusOneMatrix & rhs,
265     int numberRows, const int * whichRow,
266     int numberColumns, const int * whichColumn)
267     : ClpMatrixBase(rhs)
268{
269     matrix_ = NULL;
270     startPositive_ = NULL;
271     startNegative_ = NULL;
272     lengths_ = NULL;
273     indices_ = NULL;
274     numberRows_ = 0;
275     numberColumns_ = 0;
276     columnOrdered_ = rhs.columnOrdered_;
277#ifdef CLP_PLUS_ONE_MATRIX
278     otherFlags_ = rhs.otherFlags_;
279#endif
280     if (numberRows <= 0 || numberColumns <= 0) {
281          startPositive_ = new CoinBigIndex [1];
282          startPositive_[0] = 0;
283     } else {
284          numberColumns_ = numberColumns;
285          numberRows_ = numberRows;
286          const int * index1 = rhs.indices_;
287          CoinBigIndex * startPositive1 = rhs.startPositive_;
288
289          int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
290          int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
291          int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
292          int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
293          // Also swap incoming if not column ordered
294          if (!columnOrdered_) {
295               int temp1 = numberRows;
296               numberRows = numberColumns;
297               numberColumns = temp1;
298               const int * temp2;
299               temp2 = whichRow;
300               whichRow = whichColumn;
301               whichColumn = temp2;
302          }
303          // Throw exception if rhs empty
304          if (numberMajor1 <= 0 || numberMinor1 <= 0)
305               throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
306          // Array to say if an old row is in new copy
307          int * newRow = new int [numberMinor1];
308          int iRow;
309          for (iRow = 0; iRow < numberMinor1; iRow++)
310               newRow[iRow] = -1;
311          // and array for duplicating rows
312          int * duplicateRow = new int [numberMinor];
313          int numberBad = 0;
314          for (iRow = 0; iRow < numberMinor; iRow++) {
315               duplicateRow[iRow] = -1;
316               int kRow = whichRow[iRow];
317               if (kRow >= 0  && kRow < numberMinor1) {
318                    if (newRow[kRow] < 0) {
319                         // first time
320                         newRow[kRow] = iRow;
321                    } else {
322                         // duplicate
323                         int lastRow = newRow[kRow];
324                         newRow[kRow] = iRow;
325                         duplicateRow[iRow] = lastRow;
326                    }
327               } else {
328                    // bad row
329                    numberBad++;
330               }
331          }
332
333          if (numberBad)
334               throw CoinError("bad minor entries",
335                               "subset constructor", "ClpPlusMinusOneMatrix");
336          // now get size and check columns
337          CoinBigIndex size = 0;
338          int iColumn;
339          numberBad = 0;
340          for (iColumn = 0; iColumn < numberMajor; iColumn++) {
341               int kColumn = whichColumn[iColumn];
342               if (kColumn >= 0  && kColumn < numberMajor1) {
343                    CoinBigIndex i;
344                    for (i = startPositive1[kColumn]; i < startPositive1[kColumn+1]; i++) {
345                         int kRow = index1[i];
346                         kRow = newRow[kRow];
347                         while (kRow >= 0) {
348                              size++;
349                              kRow = duplicateRow[kRow];
350                         }
351                    }
352               } else {
353                    // bad column
354                    numberBad++;
355                    printf("%d %d %d %d\n", iColumn, numberMajor, numberMajor1, kColumn);
356               }
357          }
358          if (numberBad)
359               throw CoinError("bad major entries",
360                               "subset constructor", "ClpPlusMinusOneMatrix");
361          // now create arrays
362          startPositive_ = new CoinBigIndex [numberMajor+1];
363          startNegative_ = new CoinBigIndex [numberMajor];
364          indices_ = new int[size];
365          // and fill them
366          size = 0;
367          startPositive_[0] = 0;
368          CoinBigIndex * startNegative1 = rhs.startNegative_;
369          for (iColumn = 0; iColumn < numberMajor; iColumn++) {
370               int kColumn = whichColumn[iColumn];
371               CoinBigIndex i;
372               for (i = startPositive1[kColumn]; i < startNegative1[kColumn]; i++) {
373                    int kRow = index1[i];
374                    kRow = newRow[kRow];
375                    while (kRow >= 0) {
376                         indices_[size++] = kRow;
377                         kRow = duplicateRow[kRow];
378                    }
379               }
380               startNegative_[iColumn] = size;
381               for (; i < startPositive1[kColumn+1]; i++) {
382                    int kRow = index1[i];
383                    kRow = newRow[kRow];
384                    while (kRow >= 0) {
385                         indices_[size++] = kRow;
386                         kRow = duplicateRow[kRow];
387                    }
388               }
389               startPositive_[iColumn+1] = size;
390          }
391          delete [] newRow;
392          delete [] duplicateRow;
393     }
394     // Check valid
395     checkValid(false);
396}
397
398
399/* Returns a new matrix in reverse order without gaps */
400ClpMatrixBase *
401ClpPlusMinusOneMatrix::reverseOrderedCopy() const
402{
403     int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
404     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
405     // count number in each row/column
406     CoinBigIndex * tempP = new CoinBigIndex [numberMinor];
407     CoinBigIndex * tempN = new CoinBigIndex [numberMinor];
408     memset(tempP, 0, numberMinor * sizeof(CoinBigIndex));
409     memset(tempN, 0, numberMinor * sizeof(CoinBigIndex));
410     CoinBigIndex j = 0;
411     int i;
412     for (i = 0; i < numberMajor; i++) {
413          for (; j < startNegative_[i]; j++) {
414               int iRow = indices_[j];
415               tempP[iRow]++;
416          }
417          for (; j < startPositive_[i+1]; j++) {
418               int iRow = indices_[j];
419               tempN[iRow]++;
420          }
421     }
422     int * newIndices = new int [startPositive_[numberMajor]];
423     CoinBigIndex * newP = new CoinBigIndex [numberMinor+1];
424     CoinBigIndex * newN = new CoinBigIndex[numberMinor];
425     int iRow;
426     j = 0;
427     // do starts
428     for (iRow = 0; iRow < numberMinor; iRow++) {
429          newP[iRow] = j;
430          j += tempP[iRow];
431          tempP[iRow] = newP[iRow];
432          newN[iRow] = j;
433          j += tempN[iRow];
434          tempN[iRow] = newN[iRow];
435     }
436     newP[numberMinor] = j;
437     j = 0;
438     for (i = 0; i < numberMajor; i++) {
439          for (; j < startNegative_[i]; j++) {
440               int iRow = indices_[j];
441               CoinBigIndex put = tempP[iRow];
442               newIndices[put++] = i;
443               tempP[iRow] = put;
444          }
445          for (; j < startPositive_[i+1]; j++) {
446               int iRow = indices_[j];
447               CoinBigIndex put = tempN[iRow];
448               newIndices[put++] = i;
449               tempN[iRow] = put;
450          }
451     }
452     delete [] tempP;
453     delete [] tempN;
454     ClpPlusMinusOneMatrix * newCopy = new ClpPlusMinusOneMatrix();
455     newCopy->passInCopy(numberMinor, numberMajor,
456                         !columnOrdered_,  newIndices, newP, newN);
457     return newCopy;
458}
459//static bool doPlusOnes=true;
460//unscaled versions
461void
462ClpPlusMinusOneMatrix::times(double scalar,
463                             const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
464{
465     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
466     int i;
467     CoinBigIndex j;
468     assert (columnOrdered_);
469#ifdef CLP_PLUS_ONE_MATRIX
470     if ((otherFlags_&1)==0||!doPlusOnes) {
471#endif
472       for (i = 0; i < numberMajor; i++) {
473          double value = scalar * x[i];
474          if (value) {
475               for (j = startPositive_[i]; j < startNegative_[i]; j++) {
476                    int iRow = indices_[j];
477                    y[iRow] += value;
478               }
479               for (; j < startPositive_[i+1]; j++) {
480                    int iRow = indices_[j];
481                    y[iRow] -= value;
482               }
483          }
484       }
485#ifdef CLP_PLUS_ONE_MATRIX
486     } else {
487       // plus one
488       oneit(0);
489       for (i = 0; i < numberMajor; i++) {
490          double value = scalar * x[i];
491          if (value) {
492               for (j = startPositive_[i]; j < startPositive_[i+1]; j++) {
493                    int iRow = indices_[j];
494                    y[iRow] += value;
495               }
496          }
497       }
498     }
499#endif
500}
501void
502ClpPlusMinusOneMatrix::transposeTimes(double scalar,
503                                      const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
504{
505     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
506     int i;
507     CoinBigIndex j = 0;
508     assert (columnOrdered_);
509#ifdef CLP_PLUS_ONE_MATRIX
510     if ((otherFlags_&1)==0||!doPlusOnes) {
511#endif
512       for (i = 0; i < numberMajor; i++) {
513          double value = 0.0;
514          for (; j < startNegative_[i]; j++) {
515               int iRow = indices_[j];
516               value += x[iRow];
517          }
518          for (; j < startPositive_[i+1]; j++) {
519               int iRow = indices_[j];
520               value -= x[iRow];
521          }
522          y[i] += scalar * value;
523       }
524#ifdef CLP_PLUS_ONE_MATRIX
525     } else {
526       // plus one
527       oneit(1);
528       for (i = 0; i < numberMajor; i++) {
529          double value = 0.0;
530          for (; j < startPositive_[i+1]; j++) {
531               int iRow = indices_[j];
532               value += x[iRow];
533          }
534          y[i] += scalar * value;
535       }
536     }
537#endif
538}
539void
540ClpPlusMinusOneMatrix::times(double scalar,
541                             const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
542                             const double * /*rowScale*/,
543                             const double * /*columnScale*/) const
544{
545     // we know it is not scaled
546     times(scalar, x, y);
547}
548void
549ClpPlusMinusOneMatrix::transposeTimes( double scalar,
550                                       const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
551                                       const double * /*rowScale*/,
552                                       const double * /*columnScale*/,
553                                       double * /*spare*/) const
554{
555     // we know it is not scaled
556     transposeTimes(scalar, x, y);
557}
558/* Return <code>x * A + y</code> in <code>z</code>.
559        Squashes small elements and knows about ClpSimplex */
560void
561ClpPlusMinusOneMatrix::transposeTimes(const ClpSimplex * model, double scalar,
562                                      const CoinIndexedVector * rowArray,
563                                      CoinIndexedVector * y,
564                                      CoinIndexedVector * columnArray) const
565{
566     // we know it is not scaled
567     columnArray->clear();
568     double * COIN_RESTRICT pi = rowArray->denseVector();
569     int numberNonZero = 0;
570     int * COIN_RESTRICT index = columnArray->getIndices();
571     double * COIN_RESTRICT array = columnArray->denseVector();
572     int numberInRowArray = rowArray->getNumElements();
573     // maybe I need one in OsiSimplex
574     double zeroTolerance = model->zeroTolerance();
575     int numberRows = model->numberRows();
576     bool packed = rowArray->packedMode();
577#ifndef NO_RTTI
578     ClpPlusMinusOneMatrix* rowCopy =
579          dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
580#else
581     ClpPlusMinusOneMatrix* rowCopy =
582          static_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
583#endif
584     double factor = 0.3;
585     // We may not want to do by row if there may be cache problems
586     int numberColumns = model->numberColumns();
587     // It would be nice to find L2 cache size - for moment 512K
588     // Be slightly optimistic
589     if (numberColumns * sizeof(double) > 1000000) {
590          if (numberRows * 10 < numberColumns)
591               factor = 0.1;
592          else if (numberRows * 4 < numberColumns)
593               factor = 0.15;
594          else if (numberRows * 2 < numberColumns)
595               factor = 0.2;
596     }
597     if (numberInRowArray > factor * numberRows || !rowCopy) {
598          assert (!y->getNumElements());
599          // do by column
600          // Need to expand if packed mode
601          int iColumn;
602          CoinBigIndex j = 0;
603          assert (columnOrdered_);
604          if (packed) {
605               // need to expand pi into y
606               assert(y->capacity() >= numberRows);
607               double * COIN_RESTRICT piOld = pi;
608               pi = y->denseVector();
609               const int * COIN_RESTRICT whichRow = rowArray->getIndices();
610               int i;
611               // modify pi so can collapse to one loop
612               for (i = 0; i < numberInRowArray; i++) {
613                    int iRow = whichRow[i];
614                    pi[iRow] = scalar * piOld[i];
615               }
616#ifdef CLP_PLUS_ONE_MATRIX
617               if ((otherFlags_&1)==0||!doPlusOnes) {
618#endif
619                 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
620                    double value = 0.0;
621                    for (; j < startNegative_[iColumn]; j++) {
622                         int iRow = indices_[j];
623                         value += pi[iRow];
624                    }
625                    for (; j < startPositive_[iColumn+1]; j++) {
626                         int iRow = indices_[j];
627                         value -= pi[iRow];
628                    }
629                    if (fabs(value) > zeroTolerance) {
630                         array[numberNonZero] = value;
631                         index[numberNonZero++] = iColumn;
632                    }
633                 }
634#ifdef CLP_PLUS_ONE_MATRIX
635               } else {
636                 // plus one
637                 oneit(2);
638                 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
639                    double value = 0.0;
640                    for (; j < startPositive_[iColumn+1]; j++) {
641                         int iRow = indices_[j];
642                         value += pi[iRow];
643                    }
644                    if (fabs(value) > zeroTolerance) {
645                         array[numberNonZero] = value;
646                         index[numberNonZero++] = iColumn;
647                    }
648                 }
649               }
650#endif
651               for (i = 0; i < numberInRowArray; i++) {
652                    int iRow = whichRow[i];
653                    pi[iRow] = 0.0;
654               }
655          } else {
656                 for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
657                    double value = 0.0;
658                    for (; j < startNegative_[iColumn]; j++) {
659                         int iRow = indices_[j];
660                         value += pi[iRow];
661                    }
662                    for (; j < startPositive_[iColumn+1]; j++) {
663                         int iRow = indices_[j];
664                         value -= pi[iRow];
665                    }
666                    value *= scalar;
667                    if (fabs(value) > zeroTolerance) {
668                         index[numberNonZero++] = iColumn;
669                         array[iColumn] = value;
670                    }
671                 }
672          }
673          columnArray->setNumElements(numberNonZero);
674     } else {
675          // do by row
676          rowCopy->transposeTimesByRow(model, scalar, rowArray, y, columnArray);
677     }
678}
679/* Return <code>x * A + y</code> in <code>z</code>.
680        Squashes small elements and knows about ClpSimplex */
681void
682ClpPlusMinusOneMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
683          const CoinIndexedVector * rowArray,
684          CoinIndexedVector * y,
685          CoinIndexedVector * columnArray) const
686{
687     columnArray->clear();
688     double * COIN_RESTRICT pi = rowArray->denseVector();
689     int numberNonZero = 0;
690     int * COIN_RESTRICT index = columnArray->getIndices();
691     double * COIN_RESTRICT array = columnArray->denseVector();
692     int numberInRowArray = rowArray->getNumElements();
693     // maybe I need one in OsiSimplex
694     double zeroTolerance = model->zeroTolerance();
695     const int * COIN_RESTRICT column = indices_;
696     const CoinBigIndex * COIN_RESTRICT startPositive = startPositive_;
697     const CoinBigIndex * COIN_RESTRICT startNegative = startNegative_;
698     const int * COIN_RESTRICT whichRow = rowArray->getIndices();
699     bool packed = rowArray->packedMode();
700     if (numberInRowArray > 2) {
701          // do by rows
702          int iRow;
703          double * COIN_RESTRICT markVector = y->denseVector(); // probably empty .. but
704          int numberOriginal = 0;
705          int i;
706          if (packed) {
707               CoinBigIndex numberCovered=0;
708               int numberColumns = getNumCols();
709               bool sparse=true;
710               int target=1*numberColumns;
711               for (int i = 0; i < numberInRowArray; i++) {
712                 int iRow = whichRow[i];
713                 numberCovered += startPositive[iRow+1] - startPositive[iRow];
714                 if (numberCovered>target) {
715                   sparse=false;
716                   break;
717                 }
718               }
719               numberNonZero = 0;
720               if (sparse) {
721                 // and set up mark as char array
722                 char * COIN_RESTRICT marked = reinterpret_cast<char *> (index + columnArray->capacity());
723                 double * COIN_RESTRICT array2 = y->denseVector();
724#ifdef CLP_DEBUG
725                 int numberColumns = model->numberColumns();
726                 for (int i = 0; i < numberColumns; i++) {
727                    assert(!marked[i]);
728                    assert(!array2[i]);
729                 }
730#endif
731#ifdef CLP_PLUS_ONE_MATRIX
732                 if ((otherFlags_&1)==0||!doPlusOnes) {
733#endif
734                   for (int i = 0; i < numberInRowArray; i++) {
735                     iRow = whichRow[i];
736                     double value = pi[i] * scalar;
737                     CoinBigIndex j;
738                     for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
739                         int iColumn = column[j];
740                         if (!marked[iColumn]) {
741                              marked[iColumn] = 1;
742                              index[numberNonZero++] = iColumn;
743                         }
744                         array2[iColumn] += value;
745                      }
746                      for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
747                         int iColumn = column[j];
748                         if (!marked[iColumn]) {
749                              marked[iColumn] = 1;
750                              index[numberNonZero++] = iColumn;
751                         }
752                         array2[iColumn] -= value;
753                      }
754                    }
755#ifdef CLP_PLUS_ONE_MATRIX
756                 } else {
757                    // plus one
758                    oneit(4);
759                   for (int i = 0; i < numberInRowArray; i++) {
760                     iRow = whichRow[i];
761                     double value = pi[i] * scalar;
762                     CoinBigIndex j;
763                     for (j = startPositive[iRow]; j < startPositive[iRow+1]; j++) {
764                         int iColumn = column[j];
765                         if (!marked[iColumn]) {
766                              marked[iColumn] = 1;
767                              index[numberNonZero++] = iColumn;
768                         }
769                         array2[iColumn] += value;
770                      }
771                    }
772                }
773#endif
774                // get rid of tiny values and zero out marked
775                 numberOriginal = numberNonZero;
776                 numberNonZero = 0;
777                 for (int i = 0; i < numberOriginal; i++) {
778                    int iColumn = index[i];
779                    if (marked[iColumn]) {
780                         double value = array2[iColumn];
781                         array2[iColumn] = 0.0;
782                         marked[iColumn] = 0;
783                         if (fabs(value) > zeroTolerance) {
784                              array[numberNonZero] = value;
785                              index[numberNonZero++] = iColumn;
786                         }
787                    }
788                 }
789               } else {
790                 // not sparse
791#ifdef CLP_PLUS_ONE_MATRIX
792                 if ((otherFlags_&1)==0||!doPlusOnes) {
793#endif
794                   for (int i = 0; i < numberInRowArray; i++) {
795                     iRow = whichRow[i];
796                     double value = pi[i] * scalar;
797                     CoinBigIndex j;
798                     for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
799                         int iColumn = column[j];
800                         array[iColumn] += value;
801                     }
802                     for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
803                         int iColumn = column[j];
804                         array[iColumn] -= value;
805                     }
806                   }
807#ifdef CLP_PLUS_ONE_MATRIX
808                } else {
809                  // plus one
810                  oneit(5);
811                   for (int i = 0; i < numberInRowArray; i++) {
812                     iRow = whichRow[i];
813                     double value = pi[i] * scalar;
814                     CoinBigIndex j;
815                     for (j = startPositive[iRow]; j < startPositive[iRow+1]; j++) {
816                         int iColumn = column[j];
817                         array[iColumn] += value;
818                     }
819                   }
820                 }
821#endif
822                 // get rid of tiny values and count
823                 for (int i = 0; i < numberColumns; i++) {
824                   double value = array[i];
825                   if (value) {
826                     array[i] = 0.0;
827                     if (fabs(value) > zeroTolerance) {
828                       array[numberNonZero] = value;
829                       index[numberNonZero++] = i;
830                     }
831                   }
832                 }
833               }
834          } else {
835               numberNonZero = 0;
836               // and set up mark as char array
837               char * COIN_RESTRICT marked = reinterpret_cast<char *> (markVector);
838               for (i = 0; i < numberOriginal; i++) {
839                    int iColumn = index[i];
840                    marked[iColumn] = 0;
841               }
842               for (i = 0; i < numberInRowArray; i++) {
843                 iRow = whichRow[i];
844                 double value = pi[iRow] * scalar;
845                 CoinBigIndex j;
846                 for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
847                   int iColumn = column[j];
848                   if (!marked[iColumn]) {
849                     marked[iColumn] = 1;
850                     index[numberNonZero++] = iColumn;
851                   }
852                   array[iColumn] += value;
853                 }
854                 for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
855                   int iColumn = column[j];
856                   if (!marked[iColumn]) {
857                     marked[iColumn] = 1;
858                     index[numberNonZero++] = iColumn;
859                   }
860                   array[iColumn] -= value;
861                 }
862               }
863               // get rid of tiny values and zero out marked
864               numberOriginal = numberNonZero;
865               numberNonZero = 0;
866               for (i = 0; i < numberOriginal; i++) {
867                    int iColumn = index[i];
868                    marked[iColumn] = 0;
869                    if (fabs(array[iColumn]) > zeroTolerance) {
870                         index[numberNonZero++] = iColumn;
871                    } else {
872                         array[iColumn] = 0.0;
873                    }
874               }
875          }
876     } else if (numberInRowArray == 2) {
877          /* do by rows when two rows (do longer first when not packed
878             and shorter first if packed */
879          int iRow0 = whichRow[0];
880          int iRow1 = whichRow[1];
881          CoinBigIndex j;
882          if (packed) {
883               double pi0 = pi[0];
884               double pi1 = pi[1];
885               if (startPositive[iRow0+1] - startPositive[iRow0] >
886                         startPositive[iRow1+1] - startPositive[iRow1]) {
887                    int temp = iRow0;
888                    iRow0 = iRow1;
889                    iRow1 = temp;
890                    pi0 = pi1;
891                    pi1 = pi[0];
892               }
893               // and set up mark as char array
894               char * COIN_RESTRICT marked = reinterpret_cast<char *> (index + columnArray->capacity());
895               int * COIN_RESTRICT lookup = y->getIndices();
896               double value = pi0 * scalar;
897               int numberOriginal;
898#ifdef CLP_PLUS_ONE_MATRIX
899               if ((otherFlags_&1)==0||!doPlusOnes) {
900#endif
901                 for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
902                    int iColumn = column[j];
903                    array[numberNonZero] = value;
904                    marked[iColumn] = 1;
905                    lookup[iColumn] = numberNonZero;
906                    index[numberNonZero++] = iColumn;
907                 }
908                 for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
909                    int iColumn = column[j];
910                    array[numberNonZero] = -value;
911                    marked[iColumn] = 1;
912                    lookup[iColumn] = numberNonZero;
913                    index[numberNonZero++] = iColumn;
914                 }
915                 numberOriginal = numberNonZero;
916                 value = pi1 * scalar;
917                 for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
918                    int iColumn = column[j];
919                    if (marked[iColumn]) {
920                         int iLookup = lookup[iColumn];
921                         array[iLookup] += value;
922                    } else {
923                         if (fabs(value) > zeroTolerance) {
924                              array[numberNonZero] = value;
925                              index[numberNonZero++] = iColumn;
926                         }
927                    }
928                 }
929                 for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
930                    int iColumn = column[j];
931                    if (marked[iColumn]) {
932                         int iLookup = lookup[iColumn];
933                         array[iLookup] -= value;
934                    } else {
935                         if (fabs(value) > zeroTolerance) {
936                              array[numberNonZero] = -value;
937                              index[numberNonZero++] = iColumn;
938                         }
939                    }
940                 }
941#ifdef CLP_PLUS_ONE_MATRIX
942               } else {
943                 // plus one
944                 oneit(7);
945                 for (j = startPositive[iRow0]; j < startPositive[iRow0+1]; j++) {
946                    int iColumn = column[j];
947                    array[numberNonZero] = value;
948                    marked[iColumn] = 1;
949                    lookup[iColumn] = numberNonZero;
950                    index[numberNonZero++] = iColumn;
951                 }
952                 numberOriginal = numberNonZero;
953                 value = pi1 * scalar;
954                 for (j = startPositive[iRow1]; j < startPositive[iRow1+1]; j++) {
955                    int iColumn = column[j];
956                    if (marked[iColumn]) {
957                         int iLookup = lookup[iColumn];
958                         array[iLookup] += value;
959                    } else {
960                         if (fabs(value) > zeroTolerance) {
961                              array[numberNonZero] = value;
962                              index[numberNonZero++] = iColumn;
963                         }
964                    }
965                 }
966               }
967#endif
968               // get rid of tiny values and zero out marked
969               int nDelete = 0;
970               for (j = 0; j < numberOriginal; j++) {
971                    int iColumn = index[j];
972                    marked[iColumn] = 0;
973                    if (fabs(array[j]) <= zeroTolerance)
974                         nDelete++;
975               }
976               if (nDelete) {
977                    numberOriginal = numberNonZero;
978                    numberNonZero = 0;
979                    for (j = 0; j < numberOriginal; j++) {
980                         int iColumn = index[j];
981                         double value = array[j];
982                         array[j] = 0.0;
983                         if (fabs(value) > zeroTolerance) {
984                              array[numberNonZero] = value;
985                              index[numberNonZero++] = iColumn;
986                         }
987                    }
988               }
989          } else {
990               if (startPositive[iRow0+1] - startPositive[iRow0] <
991                         startPositive[iRow1+1] - startPositive[iRow1]) {
992                    int temp = iRow0;
993                    iRow0 = iRow1;
994                    iRow1 = temp;
995               }
996               int numberOriginal;
997               int i;
998               numberNonZero = 0;
999               double value;
1000               value = pi[iRow0] * scalar;
1001               CoinBigIndex j;
1002               for (j = startPositive[iRow0]; j < startNegative[iRow0]; j++) {
1003                 int iColumn = column[j];
1004                 index[numberNonZero++] = iColumn;
1005                 array[iColumn] = value;
1006               }
1007               for (j = startNegative[iRow0]; j < startPositive[iRow0+1]; j++) {
1008                 int iColumn = column[j];
1009                 index[numberNonZero++] = iColumn;
1010                 array[iColumn] = -value;
1011               }
1012               value = pi[iRow1] * scalar;
1013               for (j = startPositive[iRow1]; j < startNegative[iRow1]; j++) {
1014                 int iColumn = column[j];
1015                 double value2 = array[iColumn];
1016                 if (value2) {
1017                   value2 += value;
1018                 } else {
1019                   value2 = value;
1020                   index[numberNonZero++] = iColumn;
1021                 }
1022                 array[iColumn] = value2;
1023               }
1024               for (j = startNegative[iRow1]; j < startPositive[iRow1+1]; j++) {
1025                 int iColumn = column[j];
1026                 double value2 = array[iColumn];
1027                 if (value2) {
1028                   value2 -= value;
1029                 } else {
1030                   value2 = -value;
1031                   index[numberNonZero++] = iColumn;
1032                 }
1033                 array[iColumn] = value2;
1034               }
1035               // get rid of tiny values and zero out marked
1036               numberOriginal = numberNonZero;
1037               numberNonZero = 0;
1038               for (i = 0; i < numberOriginal; i++) {
1039                    int iColumn = index[i];
1040                    if (fabs(array[iColumn]) > zeroTolerance) {
1041                         index[numberNonZero++] = iColumn;
1042                    } else {
1043                         array[iColumn] = 0.0;
1044                    }
1045               }
1046          }
1047     } else if (numberInRowArray == 1) {
1048          // Just one row
1049          int iRow = rowArray->getIndices()[0];
1050          numberNonZero = 0;
1051          double value;
1052          iRow = whichRow[0];
1053          CoinBigIndex j;
1054          if (packed) {
1055               value = pi[0] * scalar;
1056               if (fabs(value) > zeroTolerance) {
1057#ifdef CLP_PLUS_ONE_MATRIX
1058                 if ((otherFlags_&1)==0||!doPlusOnes) {
1059#endif
1060                    for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
1061                         int iColumn = column[j];
1062                         array[numberNonZero] = value;
1063                         index[numberNonZero++] = iColumn;
1064                    }
1065                    for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
1066                         int iColumn = column[j];
1067                         array[numberNonZero] = -value;
1068                         index[numberNonZero++] = iColumn;
1069                    }
1070#ifdef CLP_PLUS_ONE_MATRIX
1071                 } else {
1072                   // plus one
1073                   oneit(9);
1074                    for (j = startPositive[iRow]; j < startPositive[iRow+1]; j++) {
1075                         int iColumn = column[j];
1076                         array[numberNonZero] = value;
1077                         index[numberNonZero++] = iColumn;
1078                    }
1079                 }
1080#endif
1081               }
1082          } else {
1083               value = pi[iRow] * scalar;
1084               if (fabs(value) > zeroTolerance) {
1085                  for (j = startPositive[iRow]; j < startNegative[iRow]; j++) {
1086                         int iColumn = column[j];
1087                         array[iColumn] = value;
1088                         index[numberNonZero++] = iColumn;
1089                    }
1090                    for (j = startNegative[iRow]; j < startPositive[iRow+1]; j++) {
1091                         int iColumn = column[j];
1092                         array[iColumn] = -value;
1093                         index[numberNonZero++] = iColumn;
1094                    }
1095               }
1096          }
1097     }
1098     columnArray->setNumElements(numberNonZero);
1099     if (packed)
1100          columnArray->setPacked();
1101     y->setNumElements(0);
1102}
1103/* Return <code>x *A in <code>z</code> but
1104   just for indices in y. */
1105void
1106ClpPlusMinusOneMatrix::subsetTransposeTimes(const ClpSimplex * ,
1107          const CoinIndexedVector * rowArray,
1108          const CoinIndexedVector * y,
1109          CoinIndexedVector * columnArray) const
1110{
1111     columnArray->clear();
1112     double * COIN_RESTRICT pi = rowArray->denseVector();
1113     double * COIN_RESTRICT array = columnArray->denseVector();
1114     int jColumn;
1115     int numberToDo = y->getNumElements();
1116     const int * COIN_RESTRICT which = y->getIndices();
1117     assert (!rowArray->packedMode());
1118     columnArray->setPacked();
1119#ifdef CLP_PLUS_ONE_MATRIX
1120     if ((otherFlags_&1)==0||!doPlusOnes) {
1121#endif
1122       for (jColumn = 0; jColumn < numberToDo; jColumn++) {
1123          int iColumn = which[jColumn];
1124          double value = 0.0;
1125          CoinBigIndex j = startPositive_[iColumn];
1126          for (; j < startNegative_[iColumn]; j++) {
1127               int iRow = indices_[j];
1128               value += pi[iRow];
1129          }
1130          for (; j < startPositive_[iColumn+1]; j++) {
1131               int iRow = indices_[j];
1132               value -= pi[iRow];
1133          }
1134          array[jColumn] = value;
1135       }
1136#ifdef CLP_PLUS_ONE_MATRIX
1137     } else {
1138       // plus one
1139       oneit(11);
1140       for (jColumn = 0; jColumn < numberToDo; jColumn++) {
1141          int iColumn = which[jColumn];
1142          double value = 0.0;
1143          CoinBigIndex j = startPositive_[iColumn];
1144          for (; j < startPositive_[iColumn+1]; j++) {
1145               int iRow = indices_[j];
1146               value += pi[iRow];
1147          }
1148          array[jColumn] = value;
1149       }
1150     }
1151#endif
1152}
1153/// returns number of elements in column part of basis,
1154int
1155ClpPlusMinusOneMatrix::countBasis(const int * whichColumn,
1156                                  int & numberColumnBasic)
1157{
1158     int i;
1159     CoinBigIndex numberElements = 0;
1160     for (i = 0; i < numberColumnBasic; i++) {
1161          int iColumn = whichColumn[i];
1162          numberElements += startPositive_[iColumn+1] - startPositive_[iColumn];
1163     }
1164#if COIN_BIG_INDEX
1165     if (numberElements>COIN_INT_MAX) {
1166       printf("Factorization too large\n");
1167       abort();
1168     }
1169#endif
1170     return static_cast<int>(numberElements);
1171}
1172void
1173ClpPlusMinusOneMatrix::fillBasis(ClpSimplex * ,
1174                                 const int * whichColumn,
1175                                 int & numberColumnBasic,
1176                                 int * indexRowU, int * start,
1177                                 int * rowCount, int * columnCount,
1178                                 CoinFactorizationDouble * elementU)
1179{
1180     int i;
1181     CoinBigIndex numberElements = start[0];
1182     assert (columnOrdered_);
1183     for (i = 0; i < numberColumnBasic; i++) {
1184          int iColumn = whichColumn[i];
1185          CoinBigIndex j = startPositive_[iColumn];
1186          for (; j < startNegative_[iColumn]; j++) {
1187               int iRow = indices_[j];
1188               indexRowU[numberElements] = iRow;
1189               rowCount[iRow]++;
1190               elementU[numberElements++] = 1.0;
1191          }
1192          for (; j < startPositive_[iColumn+1]; j++) {
1193               int iRow = indices_[j];
1194               indexRowU[numberElements] = iRow;
1195               rowCount[iRow]++;
1196               elementU[numberElements++] = -1.0;
1197          }
1198          start[i+1] = static_cast<int>(numberElements);
1199          columnCount[i] = static_cast<int>(numberElements - start[i]);
1200     }
1201#if COIN_BIG_INDEX
1202     if (numberElements>COIN_INT_MAX)
1203       abort();
1204#endif
1205}
1206/* Unpacks a column into an CoinIndexedvector
1207 */
1208void
1209ClpPlusMinusOneMatrix::unpack(const ClpSimplex * ,
1210                              CoinIndexedVector * rowArray,
1211                              int iColumn) const
1212{
1213     CoinBigIndex j = startPositive_[iColumn];
1214     for (; j < startNegative_[iColumn]; j++) {
1215          int iRow = indices_[j];
1216          rowArray->add(iRow, 1.0);
1217     }
1218     for (; j < startPositive_[iColumn+1]; j++) {
1219          int iRow = indices_[j];
1220          rowArray->add(iRow, -1.0);
1221     }
1222}
1223/* Unpacks a column into an CoinIndexedvector
1224** in packed foramt
1225Note that model is NOT const.  Bounds and objective could
1226be modified if doing column generation (just for this variable) */
1227void
1228ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * ,
1229                                    CoinIndexedVector * rowArray,
1230                                    int iColumn) const
1231{
1232     int * COIN_RESTRICT index = rowArray->getIndices();
1233     double * COIN_RESTRICT array = rowArray->denseVector();
1234     int number = 0;
1235     CoinBigIndex j = startPositive_[iColumn];
1236     for (; j < startNegative_[iColumn]; j++) {
1237          int iRow = indices_[j];
1238          array[number] = 1.0;
1239          index[number++] = iRow;
1240     }
1241     for (; j < startPositive_[iColumn+1]; j++) {
1242          int iRow = indices_[j];
1243          array[number] = -1.0;
1244          index[number++] = iRow;
1245     }
1246     rowArray->setNumElements(number);
1247     rowArray->setPackedMode(true);
1248}
1249/* Adds multiple of a column into an CoinIndexedvector
1250      You can use quickAdd to add to vector */
1251void
1252ClpPlusMinusOneMatrix::add(const ClpSimplex * , CoinIndexedVector * rowArray,
1253                           int iColumn, double multiplier) const
1254{
1255     CoinBigIndex j = startPositive_[iColumn];
1256     for (; j < startNegative_[iColumn]; j++) {
1257          int iRow = indices_[j];
1258          rowArray->quickAdd(iRow, multiplier);
1259     }
1260     for (; j < startPositive_[iColumn+1]; j++) {
1261          int iRow = indices_[j];
1262          rowArray->quickAdd(iRow, -multiplier);
1263     }
1264}
1265/* Adds multiple of a column into an array */
1266void
1267ClpPlusMinusOneMatrix::add(const ClpSimplex * , double * array,
1268                           int iColumn, double multiplier) const
1269{
1270     CoinBigIndex j = startPositive_[iColumn];
1271     for (; j < startNegative_[iColumn]; j++) {
1272          int iRow = indices_[j];
1273          array[iRow] += multiplier;
1274     }
1275     for (; j < startPositive_[iColumn+1]; j++) {
1276          int iRow = indices_[j];
1277          array[iRow] -= multiplier;
1278     }
1279}
1280
1281// Return a complete CoinPackedMatrix
1282CoinPackedMatrix *
1283ClpPlusMinusOneMatrix::getPackedMatrix() const
1284{
1285     if (!matrix_) {
1286          int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
1287          int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1288          CoinBigIndex numberElements = startPositive_[numberMajor];
1289          double * elements = new double [numberElements];
1290          CoinBigIndex j = 0;
1291          int i;
1292          for (i = 0; i < numberMajor; i++) {
1293               for (; j < startNegative_[i]; j++) {
1294                    elements[j] = 1.0;
1295               }
1296               for (; j < startPositive_[i+1]; j++) {
1297                    elements[j] = -1.0;
1298               }
1299          }
1300          matrix_ =  new CoinPackedMatrix(columnOrdered_, numberMinor, numberMajor,
1301                                          getNumElements(),
1302                                          elements, indices_,
1303                                          startPositive_, getVectorLengths());
1304          delete [] elements;
1305          delete [] lengths_;
1306          lengths_ = NULL;
1307     }
1308     return matrix_;
1309}
1310/* A vector containing the elements in the packed matrix. Note that there
1311   might be gaps in this list, entries that do not belong to any
1312   major-dimension vector. To get the actual elements one should look at
1313   this vector together with vectorStarts and vectorLengths. */
1314const double *
1315ClpPlusMinusOneMatrix::getElements() const
1316{
1317     if (!matrix_)
1318          getPackedMatrix();
1319     return matrix_->getElements();
1320}
1321
1322const CoinBigIndex *
1323ClpPlusMinusOneMatrix::getVectorStarts() const
1324{
1325     return startPositive_;
1326}
1327/* The lengths of the major-dimension vectors. */
1328const int *
1329ClpPlusMinusOneMatrix::getVectorLengths() const
1330{
1331     if (!lengths_) {
1332          int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1333          lengths_ = new int [numberMajor];
1334          int i;
1335          for (i = 0; i < numberMajor; i++) {
1336            lengths_[i] = static_cast<int>(startPositive_[i+1] - startPositive_[i]);
1337          }
1338     }
1339     return lengths_;
1340}
1341/* Delete the columns whose indices are listed in <code>indDel</code>. */
1342void
1343ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel)
1344{
1345     int iColumn;
1346     CoinBigIndex newSize = startPositive_[numberColumns_];;
1347     int numberBad = 0;
1348     // Use array to make sure we can have duplicates
1349     int * which = new int[numberColumns_];
1350     memset(which, 0, numberColumns_ * sizeof(int));
1351     int nDuplicate = 0;
1352     for (iColumn = 0; iColumn < numDel; iColumn++) {
1353          int jColumn = indDel[iColumn];
1354          if (jColumn < 0 || jColumn >= numberColumns_) {
1355               numberBad++;
1356          } else {
1357               newSize -= startPositive_[jColumn+1] - startPositive_[jColumn];
1358               if (which[jColumn])
1359                    nDuplicate++;
1360               else
1361                    which[jColumn] = 1;
1362          }
1363     }
1364     if (numberBad)
1365          throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
1366     int newNumber = numberColumns_ - numDel + nDuplicate;
1367     // Get rid of temporary arrays
1368     delete [] lengths_;
1369     lengths_ = NULL;
1370     delete matrix_;
1371     matrix_ = NULL;
1372     CoinBigIndex * newPositive = new CoinBigIndex [newNumber+1];
1373     CoinBigIndex * newNegative = new CoinBigIndex [newNumber];
1374     int * newIndices = new int [newSize];
1375     newNumber = 0;
1376     newSize = 0;
1377     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1378          if (!which[iColumn]) {
1379               CoinBigIndex start, end;
1380               CoinBigIndex i;
1381               start = startPositive_[iColumn];
1382               end = startNegative_[iColumn];
1383               newPositive[newNumber] = newSize;
1384               for (i = start; i < end; i++)
1385                    newIndices[newSize++] = indices_[i];
1386               start = startNegative_[iColumn];
1387               end = startPositive_[iColumn+1];
1388               newNegative[newNumber++] = newSize;
1389               for (i = start; i < end; i++)
1390                    newIndices[newSize++] = indices_[i];
1391          }
1392     }
1393     newPositive[newNumber] = newSize;
1394     delete [] which;
1395     delete [] startPositive_;
1396     startPositive_ = newPositive;
1397     delete [] startNegative_;
1398     startNegative_ = newNegative;
1399     delete [] indices_;
1400     indices_ = newIndices;
1401     numberColumns_ = newNumber;
1402}
1403/* Delete the rows whose indices are listed in <code>indDel</code>. */
1404void
1405ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel)
1406{
1407     int iRow;
1408     int numberBad = 0;
1409     // Use array to make sure we can have duplicates
1410     int * which = new int[numberRows_];
1411     memset(which, 0, numberRows_ * sizeof(int));
1412     int nDuplicate = 0;
1413     for (iRow = 0; iRow < numDel; iRow++) {
1414          int jRow = indDel[iRow];
1415          if (jRow < 0 || jRow >= numberRows_) {
1416               numberBad++;
1417          } else {
1418               if (which[jRow])
1419                    nDuplicate++;
1420               else
1421                    which[jRow] = 1;
1422          }
1423     }
1424     if (numberBad)
1425          throw CoinError("Indices out of range", "deleteRows", "ClpPlusMinusOneMatrix");
1426     CoinBigIndex iElement;
1427     CoinBigIndex numberElements = startPositive_[numberColumns_];
1428     CoinBigIndex newSize = 0;
1429     for (iElement = 0; iElement < numberElements; iElement++) {
1430          iRow = indices_[iElement];
1431          if (!which[iRow])
1432               newSize++;
1433     }
1434     int newNumber = numberRows_ - numDel + nDuplicate;
1435     // Get rid of temporary arrays
1436     delete [] lengths_;
1437     lengths_ = NULL;
1438     delete matrix_;
1439     matrix_ = NULL;
1440     // redo which
1441     int numberRows=0;
1442     for (iRow = 0; iRow < numberRows_; iRow++) {
1443       if (which[iRow]) {
1444         which[iRow]=-1; // not wanted
1445       } else {
1446         which[iRow]=numberRows;
1447         numberRows++;
1448       }
1449     }
1450     int * newIndices = new int [newSize];
1451     newSize = 0;
1452     int iColumn;
1453     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1454          CoinBigIndex start, end;
1455          CoinBigIndex i;
1456          start = startPositive_[iColumn];
1457          end = startNegative_[iColumn];
1458          startPositive_[newNumber] = newSize;
1459          for (i = start; i < end; i++) {
1460               iRow = which[indices_[i]];
1461               if (iRow>=0)
1462                    newIndices[newSize++] = iRow;
1463          }
1464          start = startNegative_[iColumn];
1465          end = startPositive_[iColumn+1];
1466          startNegative_[newNumber] = newSize;
1467          for (i = start; i < end; i++) {
1468               iRow = which[indices_[i]];
1469               if (iRow>=0)
1470                    newIndices[newSize++] = iRow;
1471          }
1472     }
1473     startPositive_[numberColumns_] = newSize;
1474     delete [] which;
1475     delete [] indices_;
1476     indices_ = newIndices;
1477     numberRows_ = newNumber;
1478}
1479bool
1480ClpPlusMinusOneMatrix::isColOrdered() const
1481{
1482     return columnOrdered_;
1483}
1484/* Number of entries in the packed matrix. */
1485CoinBigIndex
1486ClpPlusMinusOneMatrix::getNumElements() const
1487{
1488     int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
1489     if (startPositive_)
1490          return startPositive_[numberMajor];
1491     else
1492          return 0;
1493}
1494// pass in copy (object takes ownership)
1495void
1496ClpPlusMinusOneMatrix::passInCopy(int numberRows, int numberColumns,
1497                                  bool columnOrdered, int * indices,
1498                                  CoinBigIndex * startPositive, CoinBigIndex * startNegative)
1499{
1500     columnOrdered_ = columnOrdered;
1501#ifdef CLP_PLUS_ONE_MATRIX
1502     otherFlags_ = 0;
1503#endif
1504     startPositive_ = startPositive;
1505     startNegative_ = startNegative;
1506     indices_ = indices;
1507     numberRows_ = numberRows;
1508     numberColumns_ = numberColumns;
1509     // Check valid
1510     checkValid(false);
1511}
1512// Just checks matrix valid - will say if dimensions not quite right if detail
1513void
1514ClpPlusMinusOneMatrix::checkValid(bool detail) const
1515{
1516     int maxIndex = -1;
1517     int minIndex = columnOrdered_ ? numberRows_ : numberColumns_;
1518     int number = !columnOrdered_ ? numberRows_ : numberColumns_;
1519     CoinBigIndex numberElements = getNumElements();
1520     CoinBigIndex last = -1;
1521     int bad = 0;
1522     // say if all +1
1523#ifdef CLP_PLUS_ONE_MATRIX
1524     otherFlags_ = 1;
1525#endif
1526     for (int i = 0; i < number; i++) {
1527          if(startPositive_[i] < last)
1528               bad++;
1529          else
1530               last = startPositive_[i];
1531#ifdef CLP_PLUS_ONE_MATRIX
1532          if (startNegative_[i] < startPositive_[i+1])
1533            otherFlags_ &= ~1; // not all +1
1534#endif
1535          if(startNegative_[i] < last)
1536               bad++;
1537          else
1538               last = startNegative_[i];
1539     }
1540     if(startPositive_[number] < last)
1541          bad++;
1542     CoinAssertHint(!bad, "starts are not monotonic");
1543     for (CoinBigIndex cbi = 0; cbi < numberElements; cbi++) {
1544          maxIndex = CoinMax(indices_[cbi], maxIndex);
1545          minIndex = CoinMin(indices_[cbi], minIndex);
1546     }
1547     CoinAssert(maxIndex < (columnOrdered_ ? numberRows_ : numberColumns_));
1548     CoinAssert(minIndex >= 0);
1549     if (detail) {
1550          if (minIndex > 0 || maxIndex + 1 < (columnOrdered_ ? numberRows_ : numberColumns_))
1551               printf("Not full range of indices - %d to %d\n", minIndex, maxIndex);
1552     }
1553}
1554/* Given positive integer weights for each row fills in sum of weights
1555   for each column (and slack).
1556   Returns weights vector
1557*/
1558CoinBigIndex *
1559ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex * model, int * inputWeights) const
1560{
1561     int numberRows = model->numberRows();
1562     int numberColumns = model->numberColumns();
1563     int number = numberRows + numberColumns;
1564     CoinBigIndex * weights = new CoinBigIndex[number];
1565     int i;
1566     for (i = 0; i < numberColumns; i++) {
1567          CoinBigIndex j;
1568          CoinBigIndex count = 0;
1569          for (j = startPositive_[i]; j < startPositive_[i+1]; j++) {
1570               int iRow = indices_[j];
1571               count += inputWeights[iRow];
1572          }
1573          weights[i] = count;
1574     }
1575     for (i = 0; i < numberRows; i++) {
1576          weights[i+numberColumns] = inputWeights[i];
1577     }
1578     return weights;
1579}
1580// Append Columns
1581void
1582ClpPlusMinusOneMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
1583{
1584     int iColumn;
1585     CoinBigIndex size = 0;
1586     int numberBad = 0;
1587     // say if all +1
1588#ifdef CLP_PLUS_ONE_MATRIX
1589     otherFlags_ = 0;
1590#endif
1591     for (iColumn = 0; iColumn < number; iColumn++) {
1592          int n = columns[iColumn]->getNumElements();
1593          const double * element = columns[iColumn]->getElements();
1594          size += n;
1595          int i;
1596          for (i = 0; i < n; i++) {
1597               if (fabs(element[i]) != 1.0)
1598                    numberBad++;
1599          }
1600     }
1601     if (numberBad)
1602          throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix");
1603     // Get rid of temporary arrays
1604     delete [] lengths_;
1605     lengths_ = NULL;
1606     delete matrix_;
1607     matrix_ = NULL;
1608     CoinBigIndex numberNow = startPositive_[numberColumns_];
1609     CoinBigIndex * temp;
1610     temp = new CoinBigIndex [numberColumns_+1+number];
1611     CoinMemcpyN(startPositive_, (numberColumns_ + 1), temp);
1612     delete [] startPositive_;
1613     startPositive_ = temp;
1614     temp = new CoinBigIndex [numberColumns_+number];
1615     CoinMemcpyN(startNegative_, numberColumns_, temp);
1616     delete [] startNegative_;
1617     startNegative_ = temp;
1618     int * temp2 = new int [numberNow+size];
1619     CoinMemcpyN(indices_, numberNow, temp2);
1620     delete [] indices_;
1621     indices_ = temp2;
1622     // now add
1623     size = numberNow;
1624     for (iColumn = 0; iColumn < number; iColumn++) {
1625          int n = columns[iColumn]->getNumElements();
1626          const int * row = columns[iColumn]->getIndices();
1627          const double * element = columns[iColumn]->getElements();
1628          int i;
1629          for (i = 0; i < n; i++) {
1630               if (element[i] == 1.0)
1631                    indices_[size++] = row[i];
1632          }
1633          startNegative_[iColumn+numberColumns_] = size;
1634          for (i = 0; i < n; i++) {
1635               if (element[i] == -1.0)
1636                    indices_[size++] = row[i];
1637          }
1638          startPositive_[iColumn+numberColumns_+1] = size;
1639     }
1640
1641     numberColumns_ += number;
1642}
1643// Append Rows
1644void
1645ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
1646{
1647     // Allocate arrays to use for counting
1648     int * countPositive = new int [numberColumns_+1];
1649     memset(countPositive, 0, numberColumns_ * sizeof(int));
1650     int * countNegative = new int [numberColumns_];
1651     memset(countNegative, 0, numberColumns_ * sizeof(int));
1652     int iRow;
1653     CoinBigIndex size = 0;
1654     int numberBad = 0;
1655     // say if all +1
1656#ifdef CLP_PLUS_ONE_MATRIX
1657     otherFlags_ = 0;
1658#endif
1659     for (iRow = 0; iRow < number; iRow++) {
1660          int n = rows[iRow]->getNumElements();
1661          const int * column = rows[iRow]->getIndices();
1662          const double * element = rows[iRow]->getElements();
1663          size += n;
1664          int i;
1665          for (i = 0; i < n; i++) {
1666               int iColumn = column[i];
1667               if (element[i] == 1.0)
1668                    countPositive[iColumn]++;
1669               else if (element[i] == -1.0)
1670                    countNegative[iColumn]++;
1671               else
1672                    numberBad++;
1673          }
1674     }
1675     if (numberBad)
1676          throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix");
1677     // Get rid of temporary arrays
1678     delete [] lengths_;
1679     lengths_ = NULL;
1680     delete matrix_;
1681     matrix_ = NULL;
1682     CoinBigIndex numberNow = startPositive_[numberColumns_];
1683     int * newIndices = new int [numberNow+size];
1684     // Update starts and turn counts into positions
1685     // also move current indices
1686     int iColumn;
1687     CoinBigIndex numberAdded = 0;
1688     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1689          int n, move;
1690          CoinBigIndex now;
1691          now = startPositive_[iColumn];
1692          move = static_cast<int>(startNegative_[iColumn] - now);
1693          n = countPositive[iColumn];
1694          startPositive_[iColumn] += numberAdded;
1695          CoinMemcpyN(indices_+now, move, newIndices + startPositive_[iColumn]);
1696          countPositive[iColumn] = static_cast<int>(startNegative_[iColumn] + numberAdded);
1697          numberAdded += n;
1698          now = startNegative_[iColumn];
1699          move = static_cast<int>(startPositive_[iColumn+1] - now);
1700          n = countNegative[iColumn];
1701          startNegative_[iColumn] += numberAdded;
1702          CoinMemcpyN(indices_+now, move, newIndices + startNegative_[iColumn]);
1703          countNegative[iColumn] = static_cast<int>(startPositive_[iColumn+1] + numberAdded);
1704          numberAdded += n;
1705     }
1706     delete [] indices_;
1707     indices_ = newIndices;
1708     startPositive_[numberColumns_] += numberAdded;
1709     // Now put in
1710     for (iRow = 0; iRow < number; iRow++) {
1711          int newRow = numberRows_ + iRow;
1712          int n = rows[iRow]->getNumElements();
1713          const int * column = rows[iRow]->getIndices();
1714          const double * element = rows[iRow]->getElements();
1715          int i;
1716          for (i = 0; i < n; i++) {
1717               int iColumn = column[i];
1718               int put;
1719               if (element[i] == 1.0) {
1720                    put = countPositive[iColumn];
1721                    countPositive[iColumn] = put + 1;
1722               } else {
1723                    put = countNegative[iColumn];
1724                    countNegative[iColumn] = put + 1;
1725               }
1726               indices_[put] = newRow;
1727          }
1728     }
1729     delete [] countPositive;
1730     delete [] countNegative;
1731     numberRows_ += number;
1732}
1733/* Returns largest and smallest elements of both signs.
1734   Largest refers to largest absolute value.
1735*/
1736void
1737ClpPlusMinusOneMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
1738                                       double & smallestPositive, double & largestPositive)
1739{
1740     int iColumn;
1741     bool plusOne = false;
1742     bool minusOne = false;
1743     for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1744          if (startNegative_[iColumn] > startPositive_[iColumn])
1745               plusOne = true;
1746          if (startPositive_[iColumn+1] > startNegative_[iColumn])
1747               minusOne = true;
1748     }
1749     if (minusOne) {
1750          smallestNegative = -1.0;
1751          largestNegative = -1.0;
1752     } else {
1753          smallestNegative = 0.0;
1754          largestNegative = 0.0;
1755     }
1756     if (plusOne) {
1757          smallestPositive = 1.0;
1758          largestPositive = 1.0;
1759     } else {
1760          smallestPositive = 0.0;
1761          largestPositive = 0.0;
1762     }
1763}
1764// Says whether it can do partial pricing
1765bool
1766ClpPlusMinusOneMatrix::canDoPartialPricing() const
1767{
1768     return true;
1769}
1770// Partial pricing
1771void
1772ClpPlusMinusOneMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
1773                                      int & bestSequence, int & numberWanted)
1774{
1775     numberWanted = currentWanted_;
1776     int start = static_cast<int> (startFraction * numberColumns_);
1777     int end = CoinMin(static_cast<int> (endFraction * numberColumns_ + 1), numberColumns_);
1778     CoinBigIndex j;
1779     double tolerance = model->currentDualTolerance();
1780     double * COIN_RESTRICT reducedCost = model->djRegion();
1781     const double * COIN_RESTRICT duals = model->dualRowSolution();
1782     const double * COIN_RESTRICT cost = model->costRegion();
1783     double bestDj;
1784     if (bestSequence >= 0)
1785          bestDj = fabs(reducedCost[bestSequence]);
1786     else
1787          bestDj = tolerance;
1788     int sequenceOut = model->sequenceOut();
1789     int saveSequence = bestSequence;
1790     int iSequence;
1791     for (iSequence = start; iSequence < end; iSequence++) {
1792          if (iSequence != sequenceOut) {
1793               double value;
1794               ClpSimplex::Status status = model->getStatus(iSequence);
1795
1796               switch(status) {
1797
1798               case ClpSimplex::basic:
1799               case ClpSimplex::isFixed:
1800                    break;
1801               case ClpSimplex::isFree:
1802               case ClpSimplex::superBasic:
1803                    value = cost[iSequence];
1804                    j = startPositive_[iSequence];
1805                    for (; j < startNegative_[iSequence]; j++) {
1806                         int iRow = indices_[j];
1807                         value -= duals[iRow];
1808                    }
1809                    for (; j < startPositive_[iSequence+1]; j++) {
1810                         int iRow = indices_[j];
1811                         value += duals[iRow];
1812                    }
1813                    value = fabs(value);
1814                    if (value > FREE_ACCEPT * tolerance) {
1815                         numberWanted--;
1816                         // we are going to bias towards free (but only if reasonable)
1817                         value *= FREE_BIAS;
1818                         if (value > bestDj) {
1819                              // check flagged variable and correct dj
1820                              if (!model->flagged(iSequence)) {
1821                                   bestDj = value;
1822                                   bestSequence = iSequence;
1823                              } else {
1824                                   // just to make sure we don't exit before got something
1825                                   numberWanted++;
1826                              }
1827                         }
1828                    }
1829                    break;
1830               case ClpSimplex::atUpperBound:
1831                    value = cost[iSequence];
1832                    j = startPositive_[iSequence];
1833                    for (; j < startNegative_[iSequence]; j++) {
1834                         int iRow = indices_[j];
1835                         value -= duals[iRow];
1836                    }
1837                    for (; j < startPositive_[iSequence+1]; j++) {
1838                         int iRow = indices_[j];
1839                         value += duals[iRow];
1840                    }
1841                    if (value > tolerance) {
1842                         numberWanted--;
1843                         if (value > bestDj) {
1844                              // check flagged variable and correct dj
1845                              if (!model->flagged(iSequence)) {
1846                                   bestDj = value;
1847                                   bestSequence = iSequence;
1848                              } else {
1849                                   // just to make sure we don't exit before got something
1850                                   numberWanted++;
1851                              }
1852                         }
1853                    }
1854                    break;
1855               case ClpSimplex::atLowerBound:
1856                    value = cost[iSequence];
1857                    j = startPositive_[iSequence];
1858                    for (; j < startNegative_[iSequence]; j++) {
1859                         int iRow = indices_[j];
1860                         value -= duals[iRow];
1861                    }
1862                    for (; j < startPositive_[iSequence+1]; j++) {
1863                         int iRow = indices_[j];
1864                         value += duals[iRow];
1865                    }
1866                    value = -value;
1867                    if (value > tolerance) {
1868                         numberWanted--;
1869                         if (value > bestDj) {
1870                              // check flagged variable and correct dj
1871                              if (!model->flagged(iSequence)) {
1872                                   bestDj = value;
1873                                   bestSequence = iSequence;
1874                              } else {
1875                                   // just to make sure we don't exit before got something
1876                                   numberWanted++;
1877                              }
1878                         }
1879                    }
1880                    break;
1881               }
1882          }
1883          if (!numberWanted)
1884               break;
1885     }
1886     if (bestSequence != saveSequence) {
1887          // recompute dj
1888          double value = cost[bestSequence];
1889          j = startPositive_[bestSequence];
1890          for (; j < startNegative_[bestSequence]; j++) {
1891               int iRow = indices_[j];
1892               value -= duals[iRow];
1893          }
1894          for (; j < startPositive_[bestSequence+1]; j++) {
1895               int iRow = indices_[j];
1896               value += duals[iRow];
1897          }
1898          reducedCost[bestSequence] = value;
1899          savedBestSequence_ = bestSequence;
1900          savedBestDj_ = reducedCost[savedBestSequence_];
1901     }
1902     currentWanted_ = numberWanted;
1903}
1904// Allow any parts of a created CoinMatrix to be deleted
1905void
1906ClpPlusMinusOneMatrix::releasePackedMatrix() const
1907{
1908     delete matrix_;
1909     delete [] lengths_;
1910     matrix_ = NULL;
1911     lengths_ = NULL;
1912     // say if all +1
1913#ifdef CLP_PLUS_ONE_MATRIX
1914     otherFlags_ = 0;
1915#endif
1916}
1917/* Returns true if can combine transposeTimes and subsetTransposeTimes
1918   and if it would be faster */
1919bool
1920ClpPlusMinusOneMatrix::canCombine(const ClpSimplex * model,
1921                                  const CoinIndexedVector * pi) const
1922{
1923     int numberInRowArray = pi->getNumElements();
1924     int numberRows = model->numberRows();
1925     bool packed = pi->packedMode();
1926     // factor should be smaller if doing both with two pi vectors
1927     double factor = 0.27;
1928     // We may not want to do by row if there may be cache problems
1929     // It would be nice to find L2 cache size - for moment 512K
1930     // Be slightly optimistic
1931     if (numberColumns_ * sizeof(double) > 1000000) {
1932          if (numberRows * 10 < numberColumns_)
1933               factor *= 0.333333333;
1934          else if (numberRows * 4 < numberColumns_)
1935               factor *= 0.5;
1936          else if (numberRows * 2 < numberColumns_)
1937               factor *= 0.66666666667;
1938          //if (model->numberIterations()%50==0)
1939          //printf("%d nonzero\n",numberInRowArray);
1940     }
1941     // if not packed then bias a bit more towards by column
1942     if (!packed)
1943          factor *= 0.9;
1944     return (numberInRowArray > factor * numberRows || !model->rowCopy());
1945}
1946// These have to match ClpPrimalColumnSteepest version
1947#define reference(i)  (((reference[i>>5]>>(i&31))&1)!=0)
1948/* Updates two arrays for steepest and does devex weights
1949   Returns nonzero if updates reduced cost and infeas -
1950   new infeas in dj1  - This does not at present*/
1951int
1952ClpPlusMinusOneMatrix::transposeTimes2(const ClpSimplex * model,
1953                                       const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
1954                                       const CoinIndexedVector * pi2,
1955                                       CoinIndexedVector * spare,
1956                                       double * , double * ,
1957                                       double referenceIn, double devex,
1958                                       // Array for exact devex to say what is in reference framework
1959                                       unsigned int * COIN_RESTRICT reference,
1960                                       double * COIN_RESTRICT weights, double scaleFactor)
1961{
1962     // put row of tableau in dj1
1963     double * COIN_RESTRICT pi = pi1->denseVector();
1964     int numberNonZero = 0;
1965     int * COIN_RESTRICT index = dj1->getIndices();
1966     double * COIN_RESTRICT array = dj1->denseVector();
1967     int numberInRowArray = pi1->getNumElements();
1968     double zeroTolerance = model->zeroTolerance();
1969     bool packed = pi1->packedMode();
1970     // do by column
1971     int iColumn;
1972     assert (!spare->getNumElements());
1973     double * COIN_RESTRICT piWeight = pi2->denseVector();
1974     assert (!pi2->packedMode());
1975     bool killDjs = (scaleFactor == 0.0);
1976     if (!scaleFactor)
1977          scaleFactor = 1.0;
1978     // Note scale factor was -1.0
1979     if (packed) {
1980          // need to expand pi into y
1981          assert(spare->capacity() >= model->numberRows());
1982          double * COIN_RESTRICT piOld = pi;
1983          pi = spare->denseVector();
1984          const int * COIN_RESTRICT whichRow = pi1->getIndices();
1985          int i;
1986          // modify pi so can collapse to one loop
1987          for (i = 0; i < numberInRowArray; i++) {
1988               int iRow = whichRow[i];
1989               pi[iRow] = piOld[i];
1990          }
1991          CoinBigIndex j;
1992          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
1993               ClpSimplex::Status status = model->getStatus(iColumn);
1994               if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
1995               double value = 0.0;
1996               for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
1997                    int iRow = indices_[j];
1998                    value -= pi[iRow];
1999               }
2000               for (; j < startPositive_[iColumn+1]; j++) {
2001                    int iRow = indices_[j];
2002                    value += pi[iRow];
2003               }
2004               if (fabs(value) > zeroTolerance) {
2005                    // and do other array
2006                    double modification = 0.0;
2007                    for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
2008                         int iRow = indices_[j];
2009                         modification += piWeight[iRow];
2010                    }
2011                    for (; j < startPositive_[iColumn+1]; j++) {
2012                         int iRow = indices_[j];
2013                         modification -= piWeight[iRow];
2014                    }
2015                    double thisWeight = weights[iColumn];
2016                    double pivot = value * scaleFactor;
2017                    double pivotSquared = pivot * pivot;
2018                    thisWeight += pivotSquared * devex + pivot * modification;
2019                    if (thisWeight < DEVEX_TRY_NORM) {
2020                         if (referenceIn < 0.0) {
2021                              // steepest
2022                              thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2023                         } else {
2024                              // exact
2025                              thisWeight = referenceIn * pivotSquared;
2026                              if (reference(iColumn))
2027                                   thisWeight += 1.0;
2028                              thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2029                         }
2030                    }
2031                    weights[iColumn] = thisWeight;
2032                    if (!killDjs) {
2033                         array[numberNonZero] = value;
2034                         index[numberNonZero++] = iColumn;
2035                    }
2036               }
2037          }
2038          // zero out
2039          for (i = 0; i < numberInRowArray; i++) {
2040               int iRow = whichRow[i];
2041               pi[iRow] = 0.0;
2042          }
2043     } else {
2044          CoinBigIndex j;
2045          for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2046               ClpSimplex::Status status = model->getStatus(iColumn);
2047               if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
2048               double value = 0.0;
2049               for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
2050                    int iRow = indices_[j];
2051                    value -= pi[iRow];
2052               }
2053               for (; j < startPositive_[iColumn+1]; j++) {
2054                    int iRow = indices_[j];
2055                    value += pi[iRow];
2056               }
2057               if (fabs(value) > zeroTolerance) {
2058                    // and do other array
2059                    double modification = 0.0;
2060                    for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
2061                         int iRow = indices_[j];
2062                         modification += piWeight[iRow];
2063                    }
2064                    for (; j < startPositive_[iColumn+1]; j++) {
2065                         int iRow = indices_[j];
2066                         modification -= piWeight[iRow];
2067                    }
2068                    double thisWeight = weights[iColumn];
2069                    double pivot = value * scaleFactor;
2070                    double pivotSquared = pivot * pivot;
2071                    thisWeight += pivotSquared * devex + pivot * modification;
2072                    if (thisWeight < DEVEX_TRY_NORM) {
2073                         if (referenceIn < 0.0) {
2074                              // steepest
2075                              thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2076                         } else {
2077                              // exact
2078                              thisWeight = referenceIn * pivotSquared;
2079                              if (reference(iColumn))
2080                                   thisWeight += 1.0;
2081                              thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2082                         }
2083                    }
2084                    weights[iColumn] = thisWeight;
2085                    if (!killDjs) {
2086                         array[iColumn] = value;
2087                         index[numberNonZero++] = iColumn;
2088                    }
2089               }
2090          }
2091     }
2092     dj1->setNumElements(numberNonZero);
2093     spare->setNumElements(0);
2094     if (packed)
2095          dj1->setPackedMode(true);
2096     return 0;
2097}
2098// Updates second array for steepest and does devex weights
2099void
2100ClpPlusMinusOneMatrix::subsetTimes2(const ClpSimplex * ,
2101                                    CoinIndexedVector * dj1,
2102                                    const CoinIndexedVector * pi2, CoinIndexedVector *,
2103                                    double referenceIn, double devex,
2104                                    // Array for exact devex to say what is in reference framework
2105                                    unsigned int * COIN_RESTRICT reference,
2106                                    double * COIN_RESTRICT weights, double scaleFactor)
2107{
2108     int number = dj1->getNumElements();
2109     const int * COIN_RESTRICT index = dj1->getIndices();
2110     double * COIN_RESTRICT array = dj1->denseVector();
2111     assert( dj1->packedMode());
2112
2113     double * COIN_RESTRICT piWeight = pi2->denseVector();
2114     bool killDjs = (scaleFactor == 0.0);
2115     if (!scaleFactor)
2116          scaleFactor = 1.0;
2117     for (int k = 0; k < number; k++) {
2118          int iColumn = index[k];
2119          double pivot = array[k] * scaleFactor;
2120          if (killDjs)
2121               array[k] = 0.0;
2122          // and do other array
2123          double modification = 0.0;
2124          CoinBigIndex j;
2125          for (j = startPositive_[iColumn]; j < startNegative_[iColumn]; j++) {
2126               int iRow = indices_[j];
2127               modification += piWeight[iRow];
2128          }
2129          for (; j < startPositive_[iColumn+1]; j++) {
2130               int iRow = indices_[j];
2131               modification -= piWeight[iRow];
2132          }
2133          double thisWeight = weights[iColumn];
2134          double pivotSquared = pivot * pivot;
2135          thisWeight += pivotSquared * devex + pivot * modification;
2136          if (thisWeight < DEVEX_TRY_NORM) {
2137               if (referenceIn < 0.0) {
2138                    // steepest
2139                    thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
2140               } else {
2141                    // exact
2142                    thisWeight = referenceIn * pivotSquared;
2143                    if (reference(iColumn))
2144                         thisWeight += 1.0;
2145                    thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
2146               }
2147          }
2148          weights[iColumn] = thisWeight;
2149     }
2150}
2151/* Set the dimensions of the matrix. In effect, append new empty
2152   columns/rows to the matrix. A negative number for either dimension
2153   means that that dimension doesn't change. Otherwise the new dimensions
2154   MUST be at least as large as the current ones otherwise an exception
2155   is thrown. */
2156void
2157ClpPlusMinusOneMatrix::setDimensions(int newnumrows, int newnumcols)
2158{
2159     if (newnumrows < 0)
2160          newnumrows = numberRows_;
2161     if (newnumrows < numberRows_)
2162          throw CoinError("Bad new rownum (less than current)",
2163                          "setDimensions", "CoinPackedMatrix");
2164
2165     if (newnumcols < 0)
2166          newnumcols = numberColumns_;
2167     if (newnumcols < numberColumns_)
2168          throw CoinError("Bad new colnum (less than current)",
2169                          "setDimensions", "CoinPackedMatrix");
2170
2171     int number = 0;
2172     int length = 0;
2173     if (columnOrdered_) {
2174          length = numberColumns_;
2175          numberColumns_ = newnumcols;
2176          number = numberColumns_;
2177
2178     } else {
2179          length = numberRows_;
2180          numberRows_ = newnumrows;
2181          number = numberRows_;
2182     }
2183     if (number > length) {
2184          CoinBigIndex * temp;
2185          int i;
2186          CoinBigIndex end = startPositive_[length];
2187          temp = new CoinBigIndex [number+1];
2188          CoinMemcpyN(startPositive_, (length + 1), temp);
2189          delete [] startPositive_;
2190          for (i = length + 1; i < number + 1; i++)
2191               temp[i] = end;
2192          startPositive_ = temp;
2193          temp = new CoinBigIndex [number];
2194          CoinMemcpyN(startNegative_, length, temp);
2195          delete [] startNegative_;
2196          for (i = length; i < number; i++)
2197               temp[i] = end;
2198          startNegative_ = temp;
2199     }
2200}
2201#ifndef SLIM_CLP
2202/* Append a set of rows/columns to the end of the matrix. Returns number of errors
2203   i.e. if any of the new rows/columns contain an index that's larger than the
2204   number of columns-1/rows-1 (if numberOther>0) or duplicates
2205   If 0 then rows, 1 if columns */
2206int
2207ClpPlusMinusOneMatrix::appendMatrix(int number, int type,
2208                                    const CoinBigIndex * starts, const int * index,
2209                                    const double * element, int /*numberOther*/)
2210{
2211     int numberErrors = 0;
2212     // say if all +1
2213#ifdef CLP_PLUS_ONE_MATRIX
2214     otherFlags_ = 0;
2215#endif
2216     // make into CoinPackedVector
2217     CoinPackedVectorBase ** vectors =
2218          new CoinPackedVectorBase * [number];
2219     int iVector;
2220     for (iVector = 0; iVector < number; iVector++) {
2221          CoinBigIndex iStart = starts[iVector];
2222          vectors[iVector] =
2223            new CoinPackedVector(static_cast<int>(starts[iVector+1] - iStart),
2224                                    index + iStart, element + iStart);
2225     }
2226     if (type == 0) {
2227          // rows
2228          appendRows(number, vectors);
2229     } else {
2230          // columns
2231          appendCols(number, vectors);
2232     }
2233     for (iVector = 0; iVector < number; iVector++)
2234          delete vectors[iVector];
2235     delete [] vectors;
2236     return numberErrors;
2237}
2238#endif
2239#if CLP_POOL_MATRIX
2240//#############################################################################
2241// Constructors / Destructor / Assignment
2242//#############################################################################
2243
2244//-------------------------------------------------------------------
2245// Default Constructor
2246//-------------------------------------------------------------------
2247ClpPoolMatrix::ClpPoolMatrix ()
2248  : ClpMatrixBase()
2249{
2250  setType(13);
2251  matrix_ = NULL;
2252  lengths_ = NULL;
2253  elements_ = NULL;
2254  columnStart_ = NULL;
2255  stuff_ = NULL;
2256  numberRows_ = 0;
2257  numberColumns_ = 0;
2258  numberDifferent_ = 0;
2259}
2260
2261//-------------------------------------------------------------------
2262// Copy constructor
2263//-------------------------------------------------------------------
2264ClpPoolMatrix::ClpPoolMatrix (const ClpPoolMatrix & rhs)
2265  : ClpMatrixBase(rhs)
2266{
2267  setType(13);
2268  matrix_ = NULL;
2269  lengths_ = NULL;
2270  elements_ = NULL;
2271  columnStart_ = NULL;
2272  stuff_ = NULL;
2273  numberRows_ = rhs.numberRows_;
2274  numberColumns_ = rhs.numberColumns_;
2275  numberDifferent_ = rhs.numberDifferent_;
2276  if (numberColumns_) {
2277    columnStart_ = CoinCopyOfArray(rhs.columnStart_,numberColumns_+1);
2278    CoinBigIndex numberElements = columnStart_[numberColumns_];
2279    stuff_ = CoinCopyOfArray(rhs.stuff_,numberElements);
2280    elements_ = CoinCopyOfArray(rhs.elements_,numberDifferent_);
2281  }
2282}
2283static const int mmult[] = {
2284  262139, 259459, 256889, 254291, 251701, 249133, 246709, 244247};
2285inline int hashit(double value)
2286{
2287  const unsigned char * chars = reinterpret_cast<unsigned char *>(&value);
2288  int n = 0;
2289  for (int j=0;j<8;j++) 
2290    n += mmult[j]*chars[j];
2291  n = n % (1<<CLP_POOL_SIZE);
2292  return n;
2293}
2294ClpPoolMatrix::ClpPoolMatrix (const CoinPackedMatrix & rhs)
2295  : ClpMatrixBase()
2296{
2297  setType(13);
2298  matrix_ = NULL;
2299  lengths_ = NULL;
2300  numberDifferent_ = 0;
2301  assert (rhs.isColOrdered());
2302  // get matrix data pointers
2303  const int * row = rhs.getIndices();
2304  const CoinBigIndex * columnStart = rhs.getVectorStarts();
2305  const int * columnLength = rhs.getVectorLengths();
2306  const double * elementByColumn = rhs.getElements();
2307  numberColumns_ = rhs.getNumCols();
2308  numberRows_ = rhs.getNumRows();
2309  assert (numberRows_<(1<<CLP_POOL_MATRIX));
2310  elements_ = NULL;
2311  columnStart_ = new CoinBigIndex [numberColumns_+1];
2312  stuff_ = new poolInfo [rhs.getNumElements()];
2313  int maxPool = 1<<CLP_POOL_SIZE;
2314  double * tempDifferent  = new double [maxPool];
2315#define HASH 2
2316#if HASH > 1
2317  // for hashing
2318  typedef struct {
2319    int index, next;
2320  } CoinHashLink;
2321  CoinHashLink * hashThis = new CoinHashLink[4*maxPool];
2322  for ( int i = 0; i < 4*maxPool; i++ ) {
2323    hashThis[i].index = -1;
2324    hashThis[i].next = -1;
2325  }
2326#endif
2327  int numberElements=0;
2328  int hashDifferent=0;
2329  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2330    CoinBigIndex k;
2331    columnStart_[iColumn]=numberElements;
2332    for (k = columnStart[iColumn]; k < columnStart[iColumn] + columnLength[iColumn];
2333         k++) {
2334      int iRow=row[k];
2335      double value = elementByColumn[k];
2336      poolInfo stuff;
2337      stuff.row_=iRow;
2338#if HASH==1 || HASH ==3
2339      int j;
2340      for (j=0;j<numberDifferent_;j++) {
2341        if (value==tempDifferent[j])
2342          break;
2343      }
2344#endif
2345#if HASH==2 || HASH ==3
2346      int ipos = hashit(value);
2347      int j1;
2348      while ( true ) {
2349        j1 = hashThis[ipos].index;
2350        if (j1 == -1) {
2351          hashThis[ipos].index=numberDifferent_;
2352          j1=numberDifferent_;
2353          break;
2354        } else if (value==tempDifferent[j1]) {
2355          break;
2356        } else {
2357          int k = hashThis[ipos].next;
2358          if ( k == -1 ) {
2359            j1 = numberDifferent_;
2360            while ( true ) {
2361              ++hashDifferent;
2362              if ( hashThis[hashDifferent].index == -1 ) {
2363                break;
2364              }
2365            }
2366            hashThis[ipos].next = hashDifferent;
2367            hashThis[hashDifferent].index = numberDifferent_;
2368            break;
2369          } else {
2370            ipos = k;
2371            /* nothing worked - try it again */
2372          }
2373        }
2374      }
2375#if HASH ==2
2376      int j=j1;
2377#endif
2378#endif
2379#if HASH ==3
2380      assert (j==j1);
2381#endif
2382      if (j==numberDifferent_) {
2383        if (j==maxPool) 
2384          break; // too many
2385        tempDifferent[j]=value;
2386        numberDifferent_++;
2387      }
2388      stuff.pool_=j;
2389      stuff_[numberElements++]=stuff;
2390    }
2391  }
2392  columnStart_[numberColumns_]=numberElements;
2393#if HASH>1
2394  delete [] hashThis;
2395#endif
2396  elements_ = new double [numberDifferent_];
2397  memcpy(elements_,tempDifferent,numberDifferent_*sizeof(double));
2398  delete [] tempDifferent;
2399  if (numberDifferent_==maxPool) {
2400    delete [] stuff_;
2401    delete [] elements_;
2402    delete [] columnStart_;
2403    elements_ = NULL;
2404    columnStart_ = NULL;
2405    stuff_ = NULL;
2406    numberRows_ = -1;
2407    numberColumns_ = -1;
2408    numberDifferent_ = -1;
2409  }
2410}
2411
2412//-------------------------------------------------------------------
2413// Destructor
2414//-------------------------------------------------------------------
2415ClpPoolMatrix::~ClpPoolMatrix ()
2416{
2417  delete matrix_;
2418  delete [] elements_;
2419  delete [] columnStart_;
2420  delete [] lengths_;
2421  delete [] stuff_;
2422}
2423
2424//----------------------------------------------------------------
2425// Assignment operator
2426//-------------------------------------------------------------------
2427ClpPoolMatrix &
2428ClpPoolMatrix::operator=(const ClpPoolMatrix& rhs)
2429{
2430  if (this != &rhs) {
2431    ClpMatrixBase::operator=(rhs);
2432    delete matrix_;
2433    delete [] elements_;
2434    delete [] columnStart_;
2435    delete [] lengths_;
2436    delete [] stuff_;
2437    matrix_ = NULL;
2438    lengths_ = NULL;
2439    numberRows_ = rhs.numberRows_;
2440    numberColumns_ = rhs.numberColumns_;
2441    if (numberColumns_) {
2442      columnStart_ = CoinCopyOfArray(rhs.columnStart_,numberColumns_+1);
2443      CoinBigIndex numberElements = columnStart_[numberColumns_];
2444      stuff_ = CoinCopyOfArray(rhs.stuff_,numberElements);
2445      elements_ = CoinCopyOfArray(rhs.elements_,numberDifferent_);
2446    }
2447  }
2448  return *this;
2449}
2450// Constructor from arrays - handing over ownership
2451ClpPoolMatrix::ClpPoolMatrix(int numberColumns,CoinBigIndex * columnStart,
2452                             poolInfo * stuff, double * elements)
2453  : ClpMatrixBase()
2454{
2455  setType(13);
2456  matrix_ = NULL;
2457  lengths_ = NULL;
2458  numberColumns_ = numberColumns;
2459  numberDifferent_ = 0;
2460  numberRows_ = 0;
2461  columnStart_ = columnStart;
2462  stuff_ = stuff;
2463  elements_ = elements;
2464  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2465    CoinBigIndex k;
2466    for (k = columnStart_[iColumn]; k < columnStart_[iColumn+1];
2467         k++) {
2468      int iRow=stuff_[k].row_;
2469      int iPool=stuff_[k].pool_;
2470      numberDifferent_ = CoinMax(numberDifferent_,iPool);
2471      numberRows_ = CoinMax(numberRows_,iRow);
2472    }
2473  }
2474  // adjust
2475  numberDifferent_++;
2476  numberRows_++;
2477}
2478//-------------------------------------------------------------------
2479// Clone
2480//-------------------------------------------------------------------
2481ClpMatrixBase * ClpPoolMatrix::clone() const
2482{
2483  return new ClpPoolMatrix(*this);
2484}
2485/* Subset clone (without gaps).  Duplicates are allowed
2486   and order is as given */
2487ClpMatrixBase *
2488ClpPoolMatrix::subsetClone (int numberRows, const int * whichRows,
2489                            int numberColumns,
2490                            const int * whichColumns) const
2491{
2492  return new ClpPoolMatrix(*this, numberRows, whichRows,
2493                           numberColumns, whichColumns);
2494}
2495/* Subset constructor (without gaps).  Duplicates are allowed
2496   and order is as given */
2497ClpPoolMatrix::ClpPoolMatrix (
2498                              const ClpPoolMatrix & rhs,
2499                              int numberRows, const int * whichRow,
2500                              int numberColumns, const int * whichColumn)
2501  : ClpMatrixBase(rhs)
2502{
2503  setType(13);
2504  matrix_ = NULL;
2505  lengths_ = NULL;
2506  elements_ = NULL;
2507  columnStart_ = NULL;
2508  stuff_ = NULL;
2509  numberRows_ = numberRows;
2510  numberColumns_ = numberColumns;
2511  numberDifferent_ = 0;
2512  if (!numberRows_||!numberColumns_)
2513    return;
2514  numberDifferent_ = rhs.numberDifferent_;
2515  elements_ = CoinCopyOfArray(rhs.elements_,numberDifferent_);
2516  int numberRowsOther = rhs.numberRows_;
2517  int numberColumnsOther = rhs.numberColumns_;
2518  int * newRow = new int [numberRowsOther+numberRows_];
2519  int * duplicateRow = newRow+numberRowsOther;
2520  for (int iRow = 0; iRow < numberRowsOther; iRow++)
2521    newRow[iRow] = -1;
2522  // and array for duplicating rows
2523  int numberBad = 0;
2524  for (int iRow = 0; iRow < numberRows_; iRow++) {
2525    duplicateRow[iRow] = -1;
2526    int kRow = whichRow[iRow];
2527    if (kRow >= 0  && kRow < numberRowsOther) {
2528      if (newRow[kRow] < 0) {
2529        // first time
2530        newRow[kRow] = iRow;
2531      } else {
2532        // duplicate
2533        int lastRow = newRow[kRow];
2534        newRow[kRow] = iRow;
2535        duplicateRow[iRow] = lastRow;
2536      }
2537    } else {
2538      // bad row
2539      numberBad++;
2540    }
2541  }
2542 
2543  if (numberBad)
2544    throw CoinError("bad row entries",
2545                    "subset constructor", "ClpPoolMatrix");
2546  // now get size and check columns
2547  CoinBigIndex size = 0;
2548  numberBad = 0;
2549  const CoinBigIndex * starts = rhs.columnStart_;
2550  const poolInfo * stuff = rhs.stuff_;
2551  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2552    int kColumn = whichColumn[iColumn];
2553    if (kColumn >= 0  && kColumn < numberColumnsOther) {
2554      CoinBigIndex i;
2555      for (i = starts[kColumn]; i < starts[kColumn+1]; i++) {
2556        int kRow = stuff[i].row_;
2557        kRow = newRow[kRow];
2558        while (kRow >= 0) {
2559          size++;
2560          kRow = duplicateRow[kRow];
2561        }
2562      }
2563    } else {
2564      // bad column
2565      numberBad++;
2566    }
2567  }
2568  if (numberBad)
2569    throw CoinError("bad major entries",
2570                    "subset constructor", "ClpPoolMatrix");
2571  // now create arrays
2572  stuff_ = new poolInfo[size];
2573  columnStart_ = new CoinBigIndex [numberColumns_+1];
2574  size=0;
2575  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2576    int kColumn = whichColumn[iColumn];
2577    columnStart_[iColumn]=size;
2578    CoinBigIndex i;
2579    for (i = starts[kColumn]; i < starts[kColumn+1]; i++) {
2580      int kRow = stuff[i].row_;
2581      int iPool = stuff[i].pool_;
2582      kRow = newRow[kRow];
2583      while (kRow >= 0) {
2584        stuff_[size].row_=kRow;
2585        stuff_[size++].pool_=iPool;
2586        kRow = duplicateRow[kRow];
2587      }
2588    }
2589  }
2590  columnStart_[numberColumns_]=size;
2591  delete [] newRow;
2592}
2593// Create matrix_
2594ClpPackedMatrix * 
2595ClpPoolMatrix::createMatrix() const
2596{
2597  if (!matrix_) {
2598    CoinBigIndex numberElements = columnStart_[numberColumns_];
2599    double * elements = new double [numberElements];
2600    int * rows = new int [numberElements];
2601    bool lengthsExist = lengths_!=NULL;
2602    numberElements=0;
2603    if (!lengthsExist)
2604      lengths_ = new int[numberColumns_];
2605    for (int i = 0; i < numberColumns_; i++) {
2606      lengths_[i]=static_cast<int>(columnStart_[i+1]-columnStart_[i]);
2607      for (CoinBigIndex j = columnStart_[i]; j < columnStart_[i+1]; j++) {
2608        elements[numberElements] = elements_[stuff_[j].pool_];
2609        rows[numberElements++] = stuff_[j].row_;
2610      }
2611    }
2612    CoinPackedMatrix * matrix =  new CoinPackedMatrix(true, numberRows_, numberColumns_,
2613                                                      numberElements,
2614                                                      elements, rows,
2615                                                      columnStart_, lengths_);
2616    delete [] elements;
2617    delete [] rows;
2618    matrix_ = new ClpPackedMatrix(matrix);
2619  }
2620  return matrix_;
2621}
2622
2623
2624/* Returns a new matrix in reverse order without gaps */
2625ClpMatrixBase *
2626ClpPoolMatrix::reverseOrderedCopy() const
2627{
2628  // define to say no row copy
2629#define BIG_MATRIX
2630#ifndef BIG_MATRIX
2631  printf("XcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2632  createMatrix();
2633  return matrix_->reverseOrderedCopy();
2634#else
2635  return NULL;
2636#endif
2637}
2638#undef ABOCA_LITE
2639//unscaled versions
2640void
2641ClpPoolMatrix::times(double scalar,
2642                     const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
2643{
2644#if 0
2645  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2646  createMatrix()->times(scalar,x,y);
2647#else
2648  for (int iColumn = 0; iColumn < numberColumns_; iColumn++) {
2649    CoinBigIndex j;
2650    double value = x[iColumn];
2651    if (value) {
2652      CoinBigIndex start = columnStart_[iColumn];
2653      CoinBigIndex end = columnStart_[iColumn+1];
2654      value *= scalar;
2655      for (j = start; j < end; j++) {
2656        int iRow = stuff_[j].row_;
2657        y[iRow] += value * elements_[stuff_[j].pool_];
2658      }
2659    }
2660  }
2661#endif
2662}
2663void
2664ClpPoolMatrix::transposeTimes(double scalar,
2665                              const double * COIN_RESTRICT x, double * COIN_RESTRICT y) const
2666{
2667#if 0
2668  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2669  createMatrix()->transposeTimes(scalar,x,y);
2670#else
2671  int iColumn;
2672  CoinBigIndex start = columnStart_[0];
2673  if (scalar == -1.0) {
2674    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2675      CoinBigIndex j;
2676      CoinBigIndex next = columnStart_[iColumn+1];
2677      double value = 0.0;
2678      // scaled
2679      for (j = start; j < next; j++) {
2680        int jRow = stuff_[j].row_;
2681        value += x[jRow] * elements_[stuff_[j].pool_];
2682      }
2683      start = next;
2684      y[iColumn] -= value;
2685    }
2686  } else {
2687    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2688      CoinBigIndex j;
2689      CoinBigIndex next = columnStart_[iColumn+1];
2690      double value = 0.0;
2691      // scaled
2692      for (j = start; j < next; j++) {
2693        int jRow = stuff_[j].row_;
2694        value += x[jRow] * elements_[stuff_[j].pool_];
2695      }
2696      start = next;
2697      y[iColumn] += value * scalar;
2698    }
2699  }
2700#endif
2701}
2702void
2703ClpPoolMatrix::times(double scalar,
2704                     const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
2705                     const double * rowScale,
2706                     const double * columnScale) const
2707{
2708  printf("createMatrix at file %s line %d\n",__FILE__,__LINE__);
2709  createMatrix()->times(scalar,x,y);
2710}
2711void
2712ClpPoolMatrix::transposeTimes( double scalar,
2713                               const double * COIN_RESTRICT x, double * COIN_RESTRICT y,
2714                               const double * rowScale,
2715                               const double * columnScale,
2716                               double * spare) const
2717{
2718#if 0
2719  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2720  createMatrix()->transposeTimes(scalar,x,y,rowScale,columnScale,spare);
2721#else
2722  if (rowScale) {
2723    int iColumn;
2724    if (!spare) {
2725      CoinBigIndex start = columnStart_[0];
2726      if (scalar == -1.0) {
2727        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2728          CoinBigIndex j;
2729          CoinBigIndex next = columnStart_[iColumn+1];
2730          double value = 0.0;
2731          // scaled
2732          for (j = start; j < next; j++) {
2733            int jRow = stuff_[j].row_;
2734            value += x[jRow] * elements_[stuff_[j].pool_] * rowScale[jRow];
2735          }
2736          start = next;
2737          y[iColumn] -= value * columnScale[iColumn];
2738        }
2739      } else {
2740        for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2741          CoinBigIndex j;
2742          CoinBigIndex next = columnStart_[iColumn+1];
2743          double value = 0.0;
2744          // scaled
2745          for (j = start; j < next; j++) {
2746            int jRow = stuff_[j].row_;
2747            value += x[jRow] * elements_[stuff_[j].pool_] * rowScale[jRow];
2748          }
2749          start = next;
2750          y[iColumn] += value * scalar * columnScale[iColumn];
2751        }
2752      }
2753    } else {
2754      // can use spare region
2755      int iRow;
2756      int numberRows = matrix_->getNumRows();
2757      for (iRow = 0; iRow < numberRows; iRow++) {
2758        double value = x[iRow];
2759        if (value)
2760          spare[iRow] = value * rowScale[iRow];
2761        else
2762          spare[iRow] = 0.0;
2763      }
2764      CoinBigIndex start = columnStart_[0];
2765      for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
2766        CoinBigIndex j;
2767        CoinBigIndex next = columnStart_[iColumn+1];
2768        double value = 0.0;
2769        // scaled
2770        for (j = start; j < next; j++) {
2771          int jRow = stuff_[j].row_;
2772          value += spare[jRow] * elements_[stuff_[j].pool_];
2773        }
2774        start = next;
2775        y[iColumn] += value * scalar * columnScale[iColumn];
2776      }
2777    }
2778  } else {
2779    transposeTimes(scalar, x, y);
2780  }
2781#endif
2782}
2783/* Return <code>x * A + y</code> in <code>z</code>.
2784   Squashes small elements and knows about ClpSimplex */
2785void
2786ClpPoolMatrix::transposeTimes(const ClpSimplex * model, double scalar,
2787                              const CoinIndexedVector * rowArray,
2788                              CoinIndexedVector * y,
2789                              CoinIndexedVector * columnArray) const
2790{
2791#if 0
2792  bool packed = rowArray->packedMode();
2793  assert (packed);
2794  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2795  createMatrix()->transposeTimes(model,scalar,rowArray,y,columnArray);
2796#else
2797  columnArray->clear();
2798  // do by column
2799  double * COIN_RESTRICT pi = rowArray->denseVector();
2800  int numberNonZero = 0;
2801  int * COIN_RESTRICT index = columnArray->getIndices();
2802  double * COIN_RESTRICT array = columnArray->denseVector();
2803  int numberInRowArray = rowArray->getNumElements();
2804  // maybe I need one in OsiSimplex
2805  double zeroTolerance = model->zeroTolerance();
2806  bool packed = rowArray->packedMode();
2807  assert (packed);
2808  // do by column
2809  const double * COIN_RESTRICT rowScale = model->rowScale();
2810  assert (!y->getNumElements());
2811  // need to expand pi into y
2812  assert(y->capacity() >= model->numberRows());
2813  double * COIN_RESTRICT piOld = pi;
2814  pi = y->denseVector();
2815  const int * COIN_RESTRICT whichRow = rowArray->getIndices();
2816  int i;
2817  // later combine phases
2818  if (!rowScale) {
2819    // modify pi so can collapse to one loop
2820    for (i = 0; i < numberInRowArray; i++) {
2821      int iRow = whichRow[i];
2822      pi[iRow] = scalar * piOld[i];
2823    }
2824    double value = 0.0;
2825    CoinBigIndex j;
2826    CoinBigIndex end = columnStart_[1];
2827    for (j = columnStart_[0]; j < end; j++) {
2828      int iRow = stuff_[j].row_;
2829      value += pi[iRow] * elements_[stuff_[j].pool_];
2830    }
2831    int iColumn;
2832    for (iColumn = 0; iColumn < numberColumns_ - 1; iColumn++) {
2833      CoinBigIndex start = end;
2834      end = columnStart_[iColumn+2];
2835      if (fabs(value) > zeroTolerance) {
2836        array[numberNonZero] = value;
2837        index[numberNonZero++] = iColumn;
2838      }
2839      value = 0.0;
2840      for (j = start; j < end; j++) {
2841        int iRow = stuff_[j].row_;
2842        value += pi[iRow] * elements_[stuff_[j].pool_];
2843      }
2844    }
2845    if (fabs(value) > zeroTolerance) {
2846      array[numberNonZero] = value;
2847      index[numberNonZero++] = iColumn;
2848    }
2849  } else {
2850    // scaled
2851    // modify pi so can collapse to one loop
2852    for (i = 0; i < numberInRowArray; i++) {
2853      int iRow = whichRow[i];
2854      pi[iRow] = scalar * piOld[i] * rowScale[iRow];
2855    }
2856    const double * columnScale = model->columnScale();
2857    double value = 0.0;
2858    double scale = columnScale[0];
2859    CoinBigIndex j;
2860    CoinBigIndex end = columnStart_[1];
2861    for (j = columnStart_[0]; j < end; j++) {
2862      int iRow = stuff_[j].row_;
2863      value += pi[iRow] * elements_[stuff_[j].pool_];
2864    }
2865    int iColumn;
2866    for (iColumn = 0; iColumn < numberColumns_ - 1; iColumn++) {
2867      value *= scale;
2868      CoinBigIndex start = end;
2869      scale = columnScale[iColumn+1];
2870      end = columnStart_[iColumn+2];
2871      if (fabs(value) > zeroTolerance) {
2872        array[numberNonZero] = value;
2873        index[numberNonZero++] = iColumn;
2874      }
2875      value = 0.0;
2876      for (j = start; j < end; j++) {
2877        int iRow = stuff_[j].row_;
2878        value += pi[iRow] * elements_[stuff_[j].pool_];
2879      }
2880    }
2881    value *= scale;
2882    if (fabs(value) > zeroTolerance) {
2883      array[numberNonZero] = value;
2884      index[numberNonZero++] = iColumn;
2885    }
2886  }
2887  // zero out
2888  int numberRows = model->numberRows();
2889  if (numberInRowArray * 4 < numberRows) {
2890    for (i = 0; i < numberInRowArray; i++) {
2891      int iRow = whichRow[i];
2892      pi[iRow] = 0.0;
2893    }
2894  } else {
2895    CoinZeroN(pi, numberRows);
2896  }
2897  columnArray->setNumElements(numberNonZero);
2898  y->setNumElements(0);
2899  columnArray->setPackedMode(true);
2900#endif
2901}
2902/* Return <code>x * A + y</code> in <code>z</code>.
2903   Squashes small elements and knows about ClpSimplex */
2904void
2905ClpPoolMatrix::transposeTimesByRow(const ClpSimplex * model, double scalar,
2906                                   const CoinIndexedVector * rowArray,
2907                                   CoinIndexedVector * y,
2908                                   CoinIndexedVector * columnArray) const
2909{
2910  printf("createMatrix at file %s line %d\n",__FILE__,__LINE__);
2911  createMatrix()->transposeTimesByRow(model,scalar,rowArray,y,columnArray);
2912}
2913/* Return <code>x *A in <code>z</code> but
2914   just for indices in y. */
2915void
2916ClpPoolMatrix::subsetTransposeTimes(const ClpSimplex * model,
2917                                    const CoinIndexedVector * rowArray,
2918                                    const CoinIndexedVector * y,
2919                                    CoinIndexedVector * columnArray) const
2920{
2921#if 0
2922  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
2923  createMatrix()->subsetTransposeTimes(model,rowArray,y,columnArray);
2924#else
2925  columnArray->clear();
2926  double * COIN_RESTRICT pi = rowArray->denseVector();
2927  double * COIN_RESTRICT array = columnArray->denseVector();
2928  int jColumn;
2929  // get matrix data pointers
2930  const double * COIN_RESTRICT rowScale = model->rowScale();
2931  int numberToDo = y->getNumElements();
2932  const int * COIN_RESTRICT which = y->getIndices();
2933  assert (!rowArray->packedMode());
2934  columnArray->setPacked();
2935  if (!rowScale) {
2936    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
2937      int iColumn=which[jColumn];
2938      double value = 0.0;
2939      CoinBigIndex start=columnStart_[iColumn];
2940      CoinBigIndex end=columnStart_[iColumn+1];
2941      array[jColumn] = value;
2942      CoinBigIndex j;
2943      for (j = start; j < end; j++) {
2944        int iRow = stuff_[j].row_;
2945        value += pi[iRow] * elements_[stuff_[j].pool_];
2946      }
2947      array[jColumn] = value;
2948    }
2949  } else {
2950    // scaled
2951    const double * columnScale = model->columnScale();
2952    for (jColumn = 0; jColumn < numberToDo; jColumn++) {
2953      int iColumn=which[jColumn];
2954      double value = 0.0;
2955      CoinBigIndex start=columnStart_[iColumn];
2956      CoinBigIndex end=columnStart_[iColumn+1];
2957      CoinBigIndex j;
2958      for (j = start; j < end; j++) {
2959        int iRow = stuff_[j].row_;
2960        value += pi[iRow] * elements_[stuff_[j].pool_] * rowScale[iRow];
2961      }
2962      array[jColumn] = value*columnScale[iColumn];
2963    }
2964  }
2965#endif
2966}
2967/// returns number of elements in column part of basis,
2968int
2969ClpPoolMatrix::countBasis(const int * whichColumn,
2970                          int & numberColumnBasic)
2971{
2972  int i;
2973  CoinBigIndex numberElements = 0;
2974  for (i = 0; i < numberColumnBasic; i++) {
2975    int iColumn = whichColumn[i];
2976    numberElements += columnStart_[iColumn+1] - columnStart_[iColumn];
2977  }
2978#if COIN_BIG_INDEX
2979  if (numberElements>COIN_INT_MAX) {
2980    printf("Factorization too large\n");
2981    abort();
2982  }
2983#endif
2984  return static_cast<int>(numberElements);
2985}
2986void
2987ClpPoolMatrix::fillBasis(ClpSimplex * ,
2988                         const int * whichColumn,
2989                         int & numberColumnBasic,
2990                         int * indexRowU, int * start,
2991                         int * rowCount, int * columnCount,
2992                         CoinFactorizationDouble * elementU)
2993{
2994  CoinBigIndex numberElements = start[0];
2995  for (int i = 0; i < numberColumnBasic; i++) {
2996    int iColumn = whichColumn[i];
2997    CoinBigIndex j = columnStart_[iColumn];
2998    for (; j < columnStart_[iColumn+1]; j++) {
2999      int iRow = stuff_[j].row_;
3000      indexRowU[numberElements] = iRow;
3001      rowCount[iRow]++;
3002      elementU[numberElements++] = elements_[stuff_[j].pool_];;
3003    }
3004    start[i+1] = static_cast<int>(numberElements);
3005    columnCount[i] = static_cast<int>(numberElements - start[i]);
3006  }
3007#if COIN_BIG_INDEX
3008  if (numberElements>COIN_INT_MAX)
3009    abort();
3010#endif
3011}
3012/* Unpacks a column into an CoinIndexedvector
3013 */
3014void
3015ClpPoolMatrix::unpack(const ClpSimplex * ,
3016                      CoinIndexedVector * rowArray,
3017                      int iColumn) const
3018{
3019  CoinBigIndex j = columnStart_[iColumn];
3020  for (; j < columnStart_[iColumn+1]; j++) {
3021    int iRow = stuff_[j].row_;
3022    rowArray->add(iRow, elements_[stuff_[j].pool_]);
3023  }
3024}
3025/* Unpacks a column into an CoinIndexedvector
3026** in packed foramt
3027Note that model is NOT const.  Bounds and objective could
3028be modified if doing column generation (just for this variable) */
3029void
3030ClpPoolMatrix::unpackPacked(ClpSimplex * ,
3031                            CoinIndexedVector * rowArray,
3032                            int iColumn) const
3033{
3034  int * COIN_RESTRICT index = rowArray->getIndices();
3035  double * COIN_RESTRICT array = rowArray->denseVector();
3036  int number = 0;
3037  CoinBigIndex j = columnStart_[iColumn];
3038  for (; j < columnStart_[iColumn+1]; j++) {
3039    int iRow = stuff_[j].row_;
3040    array[number] = elements_[stuff_[j].pool_];
3041    index[number++] = iRow;
3042  }
3043  rowArray->setNumElements(number);
3044  rowArray->setPackedMode(true);
3045}
3046/* Adds multiple of a column into an CoinIndexedvector
3047   You can use quickAdd to add to vector */
3048void
3049ClpPoolMatrix::add(const ClpSimplex * , CoinIndexedVector * rowArray,
3050                   int iColumn, double multiplier) const
3051{
3052  CoinBigIndex j = columnStart_[iColumn];
3053  for (; j < columnStart_[iColumn+1]; j++) {
3054    int iRow = stuff_[j].row_;
3055    rowArray->quickAdd(iRow, multiplier*elements_[stuff_[j].pool_]);
3056  }
3057}
3058/* Adds multiple of a column into an array */
3059void
3060ClpPoolMatrix::add(const ClpSimplex * , double * array,
3061                   int iColumn, double multiplier) const
3062{
3063  CoinBigIndex j = columnStart_[iColumn];
3064  for (; j < columnStart_[iColumn+1]; j++) {
3065    int iRow = stuff_[j].row_;
3066    array[iRow] += multiplier*elements_[stuff_[j].pool_];
3067  }
3068}
3069
3070// Return a complete CoinPackedMatrix
3071CoinPackedMatrix *
3072ClpPoolMatrix::getPackedMatrix() const
3073{
3074  createMatrix();
3075  return matrix_->matrix();
3076}
3077/* A vector containing the elements in the packed matrix. Note that there
3078   might be gaps in this list, entries that do not belong to any
3079   major-dimension vector. To get the actual elements one should look at
3080   this vector together with vectorStarts and vectorLengths. */
3081const double *
3082ClpPoolMatrix::getElements() const
3083{
3084  printf("createMatrix at file %s line %d\n",__FILE__,__LINE__);
3085  return createMatrix()->getElements();
3086}
3087
3088/* A vector containing the minor indices of the elements in the packed
3089   matrix. Note that there might be gaps in this list, entries that do not
3090   belong to any major-dimension vector. To get the actual elements one
3091   should look at this vector together with vectorStarts and
3092   vectorLengths. */
3093const int *
3094ClpPoolMatrix::getIndices() const
3095{
3096  printf("createMatrix at file %s line %d\n",__FILE__,__LINE__);
3097  return createMatrix()->getIndices();
3098}
3099const CoinBigIndex *
3100ClpPoolMatrix::getVectorStarts() const
3101{
3102  printf("createMatrix at file %s line %d\n",__FILE__,__LINE__);
3103  return createMatrix()->getVectorStarts();
3104}
3105/* The lengths of the major-dimension vectors. */
3106const int *
3107ClpPoolMatrix::getVectorLengths() const
3108{
3109#if 0
3110  printf("XcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
3111  return createMatrix()->getVectorLengths();
3112#else
3113  if (!lengths_) {
3114    lengths_ = new int[numberColumns_];
3115    for (int i = 0; i < numberColumns_; i++)
3116      lengths_[i]=static_cast<int>(columnStart_[i+1]-columnStart_[i]);
3117  }
3118  return lengths_;
3119#endif
3120}
3121/* The length of a major-dimension vector. */
3122int
3123ClpPoolMatrix::getVectorLength(int index) const
3124{
3125  return static_cast<int>(columnStart_[index+1]-columnStart_[index]);
3126}
3127/* Delete the columns whose indices are listed in <code>indDel</code>. */
3128void
3129ClpPoolMatrix::deleteCols(const int numDel, const int * indDel)
3130{
3131  int iColumn;
3132  CoinBigIndex newSize = columnStart_[numberColumns_];;
3133  int numberBad = 0;
3134  // Use array to make sure we can have duplicates
3135  int * which = new int[numberColumns_];
3136  memset(which, 0, numberColumns_ * sizeof(int));
3137  int nDuplicate = 0;
3138  for (iColumn = 0; iColumn < numDel; iColumn++) {
3139    int jColumn = indDel[iColumn];
3140    if (jColumn < 0 || jColumn >= numberColumns_) {
3141      numberBad++;
3142    } else {
3143      newSize -= columnStart_[jColumn+1] - columnStart_[jColumn];
3144      if (which[jColumn])
3145        nDuplicate++;
3146      else
3147        which[jColumn] = 1;
3148    }
3149  }
3150  if (numberBad)
3151    throw CoinError("Indices out of range", "deleteCols", "ClpPoolMatrix");
3152  int newNumber = numberColumns_ - numDel + nDuplicate;
3153  // Get rid of temporary arrays
3154  delete matrix_;
3155  matrix_ = NULL;
3156  CoinBigIndex * columnStart = new CoinBigIndex [newNumber+1];
3157  poolInfo * stuff = new poolInfo [newSize];
3158  newNumber = 0;
3159  newSize = 0;
3160  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3161    if (!which[iColumn]) {
3162      CoinBigIndex start, end;
3163      CoinBigIndex i;
3164      start = columnStart_[iColumn];
3165      end = columnStart_[iColumn+1];
3166      columnStart[newNumber] = newSize;
3167      for (i = start; i < end; i++)
3168        stuff[newSize++] = stuff_[i];
3169    }
3170  }
3171  columnStart[newNumber] = newSize;
3172  delete [] which;
3173  delete [] columnStart_;
3174  columnStart_ = columnStart;
3175  delete [] stuff_;
3176  stuff_ = stuff;
3177  numberColumns_ = newNumber;
3178}
3179/* Delete the rows whose indices are listed in <code>indDel</code>. */
3180void
3181ClpPoolMatrix::deleteRows(const int numDel, const int * indDel)
3182{
3183  int iRow;
3184  int numberBad = 0;
3185  // Use array to make sure we can have duplicates
3186  int * which = new int[numberRows_];
3187  memset(which, 0, numberRows_ * sizeof(int));
3188  int nDuplicate = 0;
3189  for (iRow = 0; iRow < numDel; iRow++) {
3190    int jRow = indDel[iRow];
3191    if (jRow < 0 || jRow >= numberRows_) {
3192      numberBad++;
3193    } else {
3194      if (which[jRow])
3195        nDuplicate++;
3196      else
3197        which[jRow] = 1;
3198    }
3199  }
3200  if (numberBad)
3201    throw CoinError("Indices out of range", "deleteRows", "ClpPoolMatrix");
3202  CoinBigIndex iElement;
3203  CoinBigIndex numberElements = columnStart_[numberColumns_];
3204  CoinBigIndex newSize = 0;
3205  for (iElement = 0; iElement < numberElements; iElement++) {
3206    iRow = stuff_[iElement].row_;
3207    if (!which[iRow])
3208      newSize++;
3209  }
3210  int newNumber = numberRows_ - numDel + nDuplicate;
3211  // Get rid of temporary arrays
3212  delete matrix_;
3213  matrix_ = NULL;
3214  poolInfo * stuff = new poolInfo [newSize];
3215  newSize = 0;
3216  // redo which
3217  int numberRows=0;
3218  for (iRow = 0; iRow < numberRows_; iRow++) {
3219    if (which[iRow]) {
3220      which[iRow]=-1; // not wanted
3221    } else {
3222      which[iRow]=numberRows;
3223      numberRows++;
3224    }
3225  }
3226  int iColumn;
3227  for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3228    CoinBigIndex start, end;
3229    CoinBigIndex i;
3230    start = columnStart_[iColumn];
3231    end = columnStart_[iColumn+1];
3232    columnStart_[newNumber] = newSize;
3233    for (i = start; i < end; i++) {
3234      poolInfo newStuff=stuff_[i];
3235      iRow = which[newStuff.row_];
3236      if (iRow>=0) {
3237        newStuff.row_=iRow;
3238        stuff[newSize++] = newStuff;
3239      }
3240    }
3241  }
3242  columnStart_[numberColumns_] = newSize;
3243  delete [] which;
3244  delete [] stuff_;
3245  stuff_ = stuff;
3246  numberRows_ = newNumber;
3247}
3248bool
3249ClpPoolMatrix::isColOrdered() const
3250{
3251  return true;
3252}
3253/* Number of entries in the packed matrix. */
3254CoinBigIndex
3255ClpPoolMatrix::getNumElements() const
3256{
3257  if (numberColumns_)
3258    return columnStart_[numberColumns_];
3259  else
3260    return 0;
3261}
3262/* Returns largest and smallest elements of both signs.
3263   Largest refers to largest absolute value.
3264*/
3265void
3266ClpPoolMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
3267                               double & smallestPositive, double & largestPositive)
3268{
3269  smallestNegative = -COIN_DBL_MAX;
3270  largestNegative = 0.0;
3271  smallestPositive = COIN_DBL_MAX;
3272  largestPositive = 0.0;
3273  for (int i = 0; i < numberDifferent_; i++) {
3274    double value = elements_[i];
3275    if (value > 0.0) {
3276      smallestPositive = CoinMin(smallestPositive, value);
3277      largestPositive = CoinMax(largestPositive, value);
3278    } else if (value < 0.0) {
3279      smallestNegative = CoinMax(smallestNegative, value);
3280      largestNegative = CoinMin(largestNegative, value);
3281    }
3282  }
3283}
3284// Says whether it can do partial pricing
3285bool
3286ClpPoolMatrix::canDoPartialPricing() const
3287{
3288  return true;
3289}
3290// Partial pricing
3291void
3292ClpPoolMatrix::partialPricing(ClpSimplex * model, double startFraction, double endFraction,
3293                              int & bestSequence, int & numberWanted)
3294{
3295  printf("createMatrix at file %s line %d\n",__FILE__,__LINE__);
3296  createMatrix()->partialPricing(model,startFraction,
3297                                 endFraction,bestSequence,
3298                                 numberWanted);
3299}
3300// Allow any parts of a created CoinMatrix to be deleted
3301void
3302ClpPoolMatrix::releasePackedMatrix() const
3303{
3304  delete matrix_;
3305  matrix_ = NULL;
3306}
3307// These have to match ClpPrimalColumnSteepest version
3308#define reference(i)  (((reference[i>>5]>>(i&31))&1)!=0)
3309/* Updates two arrays for steepest and does devex weights
3310   Returns nonzero if updates reduced cost and infeas -
3311   new infeas in dj1  - This does not at present*/
3312int
3313ClpPoolMatrix::transposeTimes2(const ClpSimplex * model,
3314                               const CoinIndexedVector * pi1, CoinIndexedVector * dj1,
3315                               const CoinIndexedVector * pi2,
3316                               CoinIndexedVector * spare,
3317                               double * infeas, double * reducedCost,
3318                               double referenceIn, double devex,
3319                               // Array for exact devex to say what is in reference framework
3320                               unsigned int * COIN_RESTRICT reference,
3321                               double * COIN_RESTRICT weights, double scaleFactor)
3322{
3323#if 0
3324  printf("CcreateMatrix at file %s line %d\n",__FILE__,__LINE__);
3325  return createMatrix()->transposeTimes2(model,pi1,dj1,pi2,spare,
3326                                         infeas,reducedCost,referenceIn,devex,
3327                                         reference,weights,scaleFactor);
3328#else
3329  int returnCode=0;
3330  // put row of tableau in dj1
3331  double * COIN_RESTRICT pi = pi1->denseVector();
3332  int numberNonZero = 0;
3333  int * COIN_RESTRICT index = dj1->getIndices();
3334  double * COIN_RESTRICT array = dj1->denseVector();
3335  int numberInRowArray = pi1->getNumElements();
3336  double zeroTolerance = model->zeroTolerance();
3337  double dualTolerance = model->currentDualTolerance();
3338  // we can't really trust infeasibilities if there is dual error
3339  // this coding has to mimic coding in checkDualSolution
3340  double error = CoinMin(1.0e-2, model->largestDualError());
3341  // allow tolerance at least slightly bigger than standard
3342  dualTolerance = dualTolerance +  error;
3343  bool packed = pi1->packedMode();
3344  assert (packed);
3345  // do by column
3346  int iColumn;
3347  const double * COIN_RESTRICT rowScale = model->rowScale();
3348  assert (!spare->getNumElements());
3349  assert (numberColumns_ > 0);
3350  double * COIN_RESTRICT piWeight = pi2->denseVector();
3351  assert (!pi2->packedMode());
3352  bool killDjs = (scaleFactor == 0.0);
3353  if (!scaleFactor)
3354    scaleFactor = 1.0;
3355  // need to expand pi into y
3356  assert(spare->capacity() >= model->numberRows());
3357  double * COIN_RESTRICT piOld = pi;
3358  pi = spare->denseVector();
3359  const int * COIN_RESTRICT whichRow = pi1->getIndices();
3360  int i;
3361  if (!rowScale) {
3362    // modify pi so can collapse to one loop
3363    for (i = 0; i < numberInRowArray; i++) {
3364      int iRow = whichRow[i];
3365      pi[iRow] = piOld[i];
3366    }
3367    if (infeas)
3368      returnCode=1;
3369    CoinBigIndex j;
3370    CoinBigIndex end = columnStart_[0];
3371    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3372      CoinBigIndex start = end;
3373      end = columnStart_[iColumn+1];
3374      ClpSimplex::Status status = model->getStatus(iColumn);
3375      if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
3376      double value = 0.0;
3377      for (j = start; j < end; j++) {
3378        int iRow = stuff_[j].row_;
3379        value -= pi[iRow]*elements_[stuff_[j].pool_];
3380      }
3381      if (fabs(value) > zeroTolerance) {
3382        // and do other array
3383        double modification = 0.0;
3384        for (j = start; j < end; j++) {
3385          int iRow = stuff_[j].row_;
3386          modification += piWeight[iRow]*elements_[stuff_[j].pool_];
3387        }
3388        double thisWeight = weights[iColumn];
3389        double pivot = value * scaleFactor;
3390        double pivotSquared = pivot * pivot;
3391        thisWeight += pivotSquared * devex + pivot * modification;
3392        if (thisWeight < DEVEX_TRY_NORM) {
3393          if (referenceIn < 0.0) {
3394            // steepest
3395            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
3396          } else {
3397            // exact
3398            thisWeight = referenceIn * pivotSquared;
3399            if (reference(iColumn))
3400              thisWeight += 1.0;
3401            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
3402          }
3403        }
3404        weights[iColumn] = thisWeight;
3405        if (!killDjs) {
3406          value = reducedCost[iColumn]-value;
3407          reducedCost[iColumn] = value;
3408          // simplify status
3409          ClpSimplex::Status status = model->getStatus(iColumn);
3410         
3411          switch(status) {
3412           
3413          case ClpSimplex::basic:
3414          case ClpSimplex::isFixed:
3415            break;
3416          case ClpSimplex::isFree:
3417          case ClpSimplex::superBasic:
3418            if (fabs(value) > FREE_ACCEPT * dualTolerance) {
3419              // we are going to bias towards free (but only if reasonable)
3420              value *= FREE_BIAS;
3421              value *= value;
3422              // store square in list
3423              if (infeas[iColumn]) {
3424                infeas[iColumn] = value; // already there
3425              } else {
3426                array[numberNonZero]=value;
3427                index[numberNonZero++]=iColumn;
3428              }
3429            } else {
3430              array[numberNonZero]=0.0;
3431              index[numberNonZero++]=iColumn;
3432            }
3433            break;
3434          case ClpSimplex::atUpperBound:
3435            if (value > dualTolerance) {
3436              value *= value;
3437              // store square in list
3438              if (infeas[iColumn]) {
3439                infeas[iColumn] = value; // already there
3440              } else {
3441                array[numberNonZero]=value;
3442                index[numberNonZero++]=iColumn;
3443              }
3444            } else {
3445              array[numberNonZero]=0.0;
3446              index[numberNonZero++]=iColumn;
3447            }
3448            break;
3449          case ClpSimplex::atLowerBound:
3450            if (value < -dualTolerance) {
3451              value *= value;
3452              // store square in list
3453              if (infeas[iColumn]) {
3454                infeas[iColumn] = value; // already there
3455              } else {
3456                array[numberNonZero]=value;
3457                index[numberNonZero++]=iColumn;
3458              }
3459            } else {
3460              array[numberNonZero]=0.0;
3461              index[numberNonZero++]=iColumn;
3462            }
3463          }
3464        }
3465      }
3466    }
3467  } else {
3468    // scaled
3469    // modify pi so can collapse to one loop
3470    for (i = 0; i < numberInRowArray; i++) {
3471      int iRow = whichRow[i];
3472      pi[iRow] = piOld[i] * rowScale[iRow];
3473    }
3474    // can also scale piWeight as not used again
3475    int numberWeight = pi2->getNumElements();
3476    const int * indexWeight = pi2->getIndices();
3477    for (i = 0; i < numberWeight; i++) {
3478      int iRow = indexWeight[i];
3479      piWeight[iRow] *= rowScale[iRow];
3480    }
3481    if (infeas)
3482      returnCode=1;
3483    const double * COIN_RESTRICT columnScale = model->columnScale();
3484    CoinBigIndex j;
3485    CoinBigIndex end = columnStart_[0];
3486    for (iColumn = 0; iColumn < numberColumns_; iColumn++) {
3487      CoinBigIndex start = end;
3488      end = columnStart_[iColumn+1];
3489      ClpSimplex::Status status = model->getStatus(iColumn);
3490      if (status == ClpSimplex::basic || status == ClpSimplex::isFixed) continue;
3491      double scale = columnScale[iColumn];
3492      double value = 0.0;
3493      for (j = start; j < end; j++) {
3494        int iRow = stuff_[j].row_;
3495        value -= pi[iRow]*elements_[stuff_[j].pool_];
3496      }
3497      value *= scale;
3498      if (fabs(value) > zeroTolerance) {
3499        double modification = 0.0;
3500        for (j = start; j < end; j++) {
3501          int iRow = stuff_[j].row_;
3502          modification += piWeight[iRow]*elements_[stuff_[j].pool_];
3503        }
3504        modification *= scale;
3505        double thisWeight = weights[iColumn];
3506        double pivot = value * scaleFactor;
3507        double pivotSquared = pivot * pivot;
3508        thisWeight += pivotSquared * devex + pivot * modification;
3509        if (thisWeight < DEVEX_TRY_NORM) {
3510          if (referenceIn < 0.0) {
3511            // steepest
3512            thisWeight = CoinMax(DEVEX_TRY_NORM, DEVEX_ADD_ONE + pivotSquared);
3513          } else {
3514            // exact
3515            thisWeight = referenceIn * pivotSquared;
3516            if (reference(iColumn))
3517              thisWeight += 1.0;
3518            thisWeight = CoinMax(thisWeight, DEVEX_TRY_NORM);
3519          }
3520        }
3521        weights[iColumn] = thisWeight;
3522        if (!killDjs) {
3523          value = reducedCost[iColumn]-value;
3524          reducedCost[iColumn] = value;
3525          // simplify status
3526          ClpSimplex::Status status = model->getStatus(iColumn);
3527         
3528          switch(status) {
3529           
3530          case ClpSimplex::basic:
3531          case ClpSimplex::isFixed:
3532            break;
3533          case ClpSimplex::isFree:
3534          case ClpSimplex::superBasic:
3535            if (fabs(value) > FREE_ACCEPT * dualTolerance) {
3536              // we are going to bias towards free (but only if reasonable)
3537              value *= FREE_BIAS;
3538              value *= value;
3539              // store square in list
3540              if (infeas[iColumn]) {
3541                infeas[iColumn] = value; // already there
3542              } else {
3543                array[numberNonZero]=value;
3544                index[numberNonZero++]=iColumn;
3545              }
3546            } else {
3547              array[numberNonZero]=0.0;
3548              index[numberNonZero++]=iColumn;
3549            }
3550            break;
3551          case ClpSimplex::atUpperBound:
3552            if (value > dualTolerance) {
3553              value *= value;
3554              // store square in list
3555              if (infeas[iColumn]) {
3556                infeas[iColumn] = value; // already there
3557              } else {
3558                array[numberNonZero]=value;
3559                index[numberNonZero++]=iColumn;
3560              }
3561            } else {
3562              array[numberNonZero]=0.0;
3563              index[numberNonZero++]=iColumn;
3564            }
3565            break;
3566          case ClpSimplex::atLowerBound:
3567            if (value < -dualTolerance) {
3568              value *= value;
3569              // store square in list
3570              if (infeas[iColumn]) {
3571                infeas[iColumn] = value; // already there
3572              } else {
3573                array[numberNonZero]=value;
3574                index[numberNonZero++]=iColumn;
3575              }
3576            } else {
3577              array[numberNonZero]=0.0;
3578              index[numberNonZero++]=iColumn;
3579            }
3580          }
3581        }
3582      }
3583    }
3584  }
3585  // zero out
3586  int numberRows = model->numberRows();
3587  if (numberInRowArray * 4 < numberRows) {
3588    for (i = 0; i < numberInRowArray; i++) {
3589      int iRow = whichRow[i];
3590      pi[iRow] = 0.0;
3591    }
3592  } else {
3593    CoinZeroN(pi, numberRows);
3594  }
3595  dj1->setNumElements(numberNonZero);
3596  spare->setNumElements(0);
3597  if (packed)
3598    dj1->setPackedMode(true);
3599  return returnCode;
3600#endif
3601}
3602// Updates second array for steepest and does devex weights
3603void
3604ClpPoolMatrix::subsetTimes2(const ClpSimplex * model,
3605                            CoinIndexedVector * dj1,
3606                            const CoinIndexedVector * pi2, CoinIndexedVector * dj2,
3607                            double referenceIn, double devex,
3608                            // Array for exact devex to say what is in reference framework
3609                            unsigned int * COIN_RESTRICT reference,
3610                            double * COIN_RESTRICT weights, double scaleFactor)
3611{
3612  printf("createMatrix at file %s line %d\n",__FILE__,__LINE__);
3613  createMatrix()->subsetTimes2(model,dj1,pi2,dj2,
3614                               referenceIn,devex,
3615                               reference,weights,scaleFactor);
3616}
3617/* Set the dimensions of the matrix. In effect, append new empty
3618   columns/rows to the matrix. A negative number for either dimension
3619   means that that dimension doesn't change. Otherwise the new dimensions
3620   MUST be at least as large as the current ones otherwise an exception
3621   is thrown. */
3622void
3623ClpPoolMatrix::setDimensions(int newnumrows, int newnumcols)
3624{
3625  // at present abort
3626  abort();
3627}
3628#endif
Note: See TracBrowser for help on using the repository browser.