Ignore:
Timestamp:
Oct 15, 2003 5:34:57 PM (16 years ago)
Author:
forrest
Message:

This should break everything

File:
1 edited

Legend:

Unmodified
Added
Removed
  • trunk/ClpPlusMinusOneMatrix.cpp

    r124 r225  
    5151  columnOrdered_=rhs.columnOrdered_;
    5252  if (numberColumns_) {
    53     indices_ = new int [ 2*numberColumns_];
    54     memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
     53    int numberElements = rhs.startPositive_[numberColumns_];
     54    indices_ = new int [ numberElements];
     55    memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
    5556    startPositive_ = new int [ numberColumns_+1];
    5657    memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
     
    118119    delete [] indices_;
    119120    // put in message
    120     printf("Not all +-1 - can test if indices_ null\n");
    121121    indices_=NULL;
    122122    numberRows_=0;
     
    165165    columnOrdered_=rhs.columnOrdered_;
    166166    if (numberColumns_) {
    167       indices_ = new int [ 2*numberColumns_];
    168       memcpy(indices_,rhs.indices_,2*numberColumns_*sizeof(int));
     167      int numberElements = rhs.startPositive_[numberColumns_];
     168      indices_ = new int [ numberElements];
     169      memcpy(indices_,rhs.indices_,numberElements*sizeof(int));
    169170      startPositive_ = new int [ numberColumns_+1];
    170171      memcpy(startPositive_,rhs.startPositive_,(numberColumns_+1)*sizeof(int));
     
    182183  return new ClpPlusMinusOneMatrix(*this);
    183184}
     185/* Subset clone (without gaps).  Duplicates are allowed
     186   and order is as given */
     187ClpMatrixBase *
     188ClpPlusMinusOneMatrix::subsetClone (int numberRows, const int * whichRows,
     189                              int numberColumns,
     190                              const int * whichColumns) const
     191{
     192  return new ClpPlusMinusOneMatrix(*this, numberRows, whichRows,
     193                                   numberColumns, whichColumns);
     194}
     195/* Subset constructor (without gaps).  Duplicates are allowed
     196   and order is as given */
     197ClpPlusMinusOneMatrix::ClpPlusMinusOneMatrix (
     198                       const ClpPlusMinusOneMatrix & rhs,
     199                       int numberRows, const int * whichRow,
     200                       int numberColumns, const int * whichColumn)
     201: ClpMatrixBase(rhs)
     202
     203  elements_ = NULL;
     204  startPositive_ = NULL;
     205  startNegative_ = NULL;
     206  lengths_=NULL;
     207  indices_=NULL;
     208  numberRows_=0;
     209  numberColumns_=0;
     210  columnOrdered_=rhs.columnOrdered_;
     211  if (numberRows<=0||numberColumns<=0) {
     212    startPositive_ = new int[1];
     213    startPositive_[0] = 0;
     214  } else {
     215    numberColumns_ = numberColumns;
     216    numberRows_ = numberRows;
     217    const int * index1 = rhs.indices_;
     218    int * startPositive1 = rhs.startPositive_;
     219
     220    int numberMinor = (!columnOrdered_) ? numberColumns_ : numberRows_;
     221    int numberMajor = (columnOrdered_) ? numberColumns_ : numberRows_;
     222    int numberMinor1 = (!columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
     223    int numberMajor1 = (columnOrdered_) ? rhs.numberColumns_ : rhs.numberRows_;
     224    // Also swap incoming if not column ordered
     225    if (!columnOrdered_) {
     226      int temp1 = numberRows;
     227      numberRows = numberColumns;
     228      numberColumns = temp1;
     229      const int * temp2;
     230      temp2 = whichRow;
     231      whichRow = whichColumn;
     232      whichColumn = temp2;
     233    }
     234    // Throw exception if rhs empty
     235    if (numberMajor1 <= 0 || numberMinor1 <= 0)
     236      throw CoinError("empty rhs", "subset constructor", "ClpPlusMinusOneMatrix");
     237    // Array to say if an old row is in new copy
     238    int * newRow = new int [numberMinor1];
     239    int iRow;
     240    for (iRow=0;iRow<numberMinor1;iRow++)
     241      newRow[iRow] = -1;
     242    // and array for duplicating rows
     243    int * duplicateRow = new int [numberMinor];
     244    int numberBad=0;
     245    for (iRow=0;iRow<numberMinor;iRow++) {
     246      duplicateRow[iRow] = -1;
     247      int kRow = whichRow[iRow];
     248      if (kRow>=0  && kRow < numberMinor1) {
     249        if (newRow[kRow]<0) {
     250          // first time
     251          newRow[kRow]=iRow;
     252        } else {
     253          // duplicate
     254          int lastRow = newRow[kRow];
     255          newRow[kRow]=iRow;
     256          duplicateRow[iRow] = lastRow;
     257        }
     258      } else {
     259        // bad row
     260        numberBad++;
     261      }
     262    }
     263
     264    if (numberBad)
     265      throw CoinError("bad minor entries",
     266                      "subset constructor", "ClpPlusMinusOneMatrix");
     267    // now get size and check columns
     268    int size = 0;
     269    int iColumn;
     270    numberBad=0;
     271    for (iColumn=0;iColumn<numberMajor;iColumn++) {
     272      int kColumn = whichColumn[iColumn];
     273      if (kColumn>=0  && kColumn <numberMajor1) {
     274        int i;
     275        for (i=startPositive1[kColumn];i<startPositive1[kColumn+1];i++) {
     276          int kRow = index1[i];
     277          kRow = newRow[kRow];
     278          while (kRow>=0) {
     279            size++;
     280            kRow = duplicateRow[kRow];
     281          }
     282        }
     283      } else {
     284        // bad column
     285        numberBad++;
     286        printf("%d %d %d %d\n",iColumn,numberMajor,numberMajor1,kColumn);
     287      }
     288    }
     289    if (numberBad)
     290      throw CoinError("bad major entries",
     291                      "subset constructor", "ClpPlusMinusOneMatrix");
     292    // now create arrays
     293    startPositive_ = new int [numberMajor+1];
     294    startNegative_ = new int [numberMajor];
     295    indices_ = new int[size];
     296    // and fill them
     297    size = 0;
     298    startPositive_[0]=0;
     299    int * startNegative1 = rhs.startNegative_;
     300    for (iColumn=0;iColumn<numberMajor;iColumn++) {
     301      int kColumn = whichColumn[iColumn];
     302      int i;
     303      for (i=startPositive1[kColumn];i<startNegative1[kColumn];i++) {
     304        int kRow = index1[i];
     305        kRow = newRow[kRow];
     306        while (kRow>=0) {
     307          indices_[size++] = kRow;
     308          kRow = duplicateRow[kRow];
     309        }
     310      }
     311      startNegative_[iColumn] = size;
     312      for (;i<startPositive1[kColumn+1];i++) {
     313        int kRow = index1[i];
     314        kRow = newRow[kRow];
     315        while (kRow>=0) {
     316          indices_[size++] = kRow;
     317          kRow = duplicateRow[kRow];
     318        }
     319      }
     320      startPositive_[iColumn+1] = size;
     321    }
     322    delete [] newRow;
     323    delete [] duplicateRow;
     324  }
     325}
     326
    184327
    185328/* Returns a new matrix in reverse order without gaps */
     
    323466  double zeroTolerance = model->factorization()->zeroTolerance();
    324467  int numberRows = model->numberRows();
     468  bool packed = rowArray->packedMode();
    325469  ClpPlusMinusOneMatrix* rowCopy =
    326470    dynamic_cast< ClpPlusMinusOneMatrix*>(model->rowCopy());
    327   if (numberInRowArray>0.3*numberRows||!rowCopy) {
     471  double factor = 0.3;
     472  // We may not want to do by row if there may be cache problems
     473  int numberColumns = model->numberColumns();
     474  // It would be nice to find L2 cache size - for moment 512K
     475  // Be slightly optimistic
     476  if (numberColumns*sizeof(double)>1000000) {
     477    if (numberRows*10<numberColumns)
     478      factor=0.1;
     479    else if (numberRows*4<numberColumns)
     480      factor=0.15;
     481    else if (numberRows*2<numberColumns)
     482      factor=0.2;
     483  }
     484  if (numberInRowArray>factor*numberRows||!rowCopy) {
     485    assert (!y->getNumElements());
    328486    // do by column
     487    // Need to expand if packed mode
    329488    int iColumn;
    330     double * markVector = y->denseVector(); // probably empty
    331489    CoinBigIndex j=0;
    332490    assert (columnOrdered_);
    333     for (iColumn=0;iColumn<numberColumns_;iColumn++) {
    334       double value2 = 0.0;
    335       for (;j<startNegative_[iColumn];j++) {
    336         int iRow = indices_[j];
    337         value2 += pi[iRow];
    338       }
    339       for (;j<startPositive_[iColumn+1];j++) {
    340         int iRow = indices_[j];
    341         value2 -= pi[iRow];
    342       }
    343       double value = markVector[iColumn];
    344       markVector[iColumn]=0.0;
    345       value += scalar*value2;
    346       if (fabs(value)>zeroTolerance) {
    347         index[numberNonZero++]=iColumn;
    348         array[iColumn]=value;
     491    if (packed) {
     492      // need to expand pi into y
     493      assert(y->capacity()>=numberRows);
     494      double * piOld = pi;
     495      pi = y->denseVector();
     496      const int * whichRow = rowArray->getIndices();
     497      int i;
     498      // modify pi so can collapse to one loop
     499      for (i=0;i<numberInRowArray;i++) {
     500        int iRow = whichRow[i];
     501        pi[iRow]=scalar*piOld[i];
     502      }
     503      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     504        double value = 0.0;
     505        for (;j<startNegative_[iColumn];j++) {
     506          int iRow = indices_[j];
     507          value += pi[iRow];
     508        }
     509        for (;j<startPositive_[iColumn+1];j++) {
     510          int iRow = indices_[j];
     511          value -= pi[iRow];
     512        }
     513        if (fabs(value)>zeroTolerance) {
     514          array[numberNonZero]=value;
     515          index[numberNonZero++]=iColumn;
     516        }
     517      }
     518      for (i=0;i<numberInRowArray;i++) {
     519        int iRow = whichRow[i];
     520        pi[iRow]=0.0;
     521      }
     522    } else {
     523      for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     524        double value = 0.0;
     525        for (;j<startNegative_[iColumn];j++) {
     526          int iRow = indices_[j];
     527          value += pi[iRow];
     528        }
     529        for (;j<startPositive_[iColumn+1];j++) {
     530          int iRow = indices_[j];
     531          value -= pi[iRow];
     532        }
     533        value *= scalar;
     534        if (fabs(value)>zeroTolerance) {
     535          index[numberNonZero++]=iColumn;
     536          array[iColumn]=value;
     537        }
    349538      }
    350539    }
    351540    columnArray->setNumElements(numberNonZero);
    352     y->setNumElements(0);
    353541  } else {
    354542    // do by row
     
    376564  const CoinBigIndex * startNegative = startNegative_;
    377565  const int * whichRow = rowArray->getIndices();
     566  bool packed = rowArray->packedMode();
    378567  if (numberInRowArray>2||y->getNumElements()) {
    379568    // do by rows
     
    383572    int numberOriginal=y->getNumElements();
    384573    int i;
    385     for (i=0;i<numberOriginal;i++) {
    386       int iColumn = mark[i];
    387       index[i]=iColumn;
    388       array[iColumn]=markVector[iColumn];
    389       markVector[iColumn]=0.0;
    390     }
    391     numberNonZero=numberOriginal;
    392     // and set up mark as char array
    393     char * marked = (char *) markVector;
    394     for (i=0;i<numberOriginal;i++) {
    395       int iColumn = index[i];
    396       marked[iColumn]=0;
    397     }
    398     for (i=0;i<numberInRowArray;i++) {
    399       iRow = whichRow[i];
    400       double value = pi[iRow]*scalar;
    401       CoinBigIndex j;
    402       for (j=startPositive[iRow];j<startNegative[iRow];j++) {
    403         int iColumn = column[j];
    404         if (!marked[iColumn]) {
    405           marked[iColumn]=1;
     574    if (packed) {
     575      assert(!numberOriginal);
     576      numberNonZero=0;
     577      // and set up mark as char array
     578      char * marked = (char *) (index+columnArray->capacity());
     579      double * array2 = y->denseVector();
     580#ifdef CLP_DEBUG
     581      int numberColumns = model->numberColumns();
     582      for (i=0;i<numberColumns;i++) {
     583        assert(!marked[i]);
     584        assert(!array2[i]);
     585      }
     586#endif
     587      for (i=0;i<numberInRowArray;i++) {
     588        iRow = whichRow[i];
     589        double value = pi[i]*scalar;
     590        CoinBigIndex j;
     591        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     592          int iColumn = column[j];
     593          if (!marked[iColumn]) {
     594            marked[iColumn]=1;
     595            index[numberNonZero++]=iColumn;
     596          }
     597          array2[iColumn] += value;
     598        }
     599        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     600          int iColumn = column[j];
     601          if (!marked[iColumn]) {
     602            marked[iColumn]=1;
     603            index[numberNonZero++]=iColumn;
     604          }
     605          array2[iColumn] -= value;
     606        }
     607      }
     608      // get rid of tiny values and zero out marked
     609      numberOriginal=numberNonZero;
     610      numberNonZero=0;
     611      for (i=0;i<numberOriginal;i++) {
     612        int iColumn = index[i];
     613        if (marked[iColumn]) {
     614          double value = array2[iColumn];
     615          array2[iColumn]=0.0;
     616          marked[iColumn]=0;
     617          if (fabs(value)>zeroTolerance) {
     618            array[numberNonZero]=value;
     619            index[numberNonZero++]=iColumn;
     620          }
     621        }
     622      }
     623    } else {     
     624      for (i=0;i<numberOriginal;i++) {
     625        int iColumn = mark[i];
     626        index[i]=iColumn;
     627        array[iColumn]=markVector[iColumn];
     628        markVector[iColumn]=0.0;
     629      }
     630      numberNonZero=numberOriginal;
     631      // and set up mark as char array
     632      char * marked = (char *) markVector;
     633      for (i=0;i<numberOriginal;i++) {
     634        int iColumn = index[i];
     635        marked[iColumn]=0;
     636      }
     637      for (i=0;i<numberInRowArray;i++) {
     638        iRow = whichRow[i];
     639        double value = pi[iRow]*scalar;
     640        CoinBigIndex j;
     641        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     642          int iColumn = column[j];
     643          if (!marked[iColumn]) {
     644            marked[iColumn]=1;
     645            index[numberNonZero++]=iColumn;
     646          }
     647          array[iColumn] += value;
     648        }
     649        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     650          int iColumn = column[j];
     651          if (!marked[iColumn]) {
     652            marked[iColumn]=1;
     653            index[numberNonZero++]=iColumn;
     654          }
     655          array[iColumn] -= value;
     656        }
     657      }
     658      // get rid of tiny values and zero out marked
     659      numberOriginal=numberNonZero;
     660      numberNonZero=0;
     661      for (i=0;i<numberOriginal;i++) {
     662        int iColumn = index[i];
     663        marked[iColumn]=0;
     664        if (fabs(array[iColumn])>zeroTolerance) {
    406665          index[numberNonZero++]=iColumn;
    407         }
    408         array[iColumn] += value;
    409       }
    410       for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
    411         int iColumn = column[j];
    412         if (!marked[iColumn]) {
    413           marked[iColumn]=1;
    414           index[numberNonZero++]=iColumn;
    415         }
    416         array[iColumn] -= value;
    417       }
    418     }
    419     // get rid of tiny values and zero out marked
    420     numberOriginal=numberNonZero;
    421     numberNonZero=0;
    422     for (i=0;i<numberOriginal;i++) {
    423       int iColumn = index[i];
    424       marked[iColumn]=0;
    425       if (fabs(array[iColumn])>zeroTolerance) {
    426         index[numberNonZero++]=iColumn;
    427       } else {
    428         array[iColumn]=0.0;
     666        } else {
     667          array[iColumn]=0.0;
     668        }
    429669      }
    430670    }
    431671  } else if (numberInRowArray==2) {
    432     // do by rows when two rows (do longer first)
     672    /* do by rows when two rows (do longer first when not packed
     673       and shorter first if packed */
    433674    int iRow0 = whichRow[0];
    434675    int iRow1 = whichRow[1];
    435     if (startPositive[iRow0+1]-startPositive[iRow0]<
    436         startPositive[iRow1+1]-startPositive[iRow1]) {
    437       int temp = iRow0;
    438       iRow0 = iRow1;
    439       iRow1 = temp;
    440     }
    441     int numberOriginal;
    442     int i;
    443     numberNonZero=0;
    444     double value;
    445     value = pi[iRow0]*scalar;
    446     CoinBigIndex j;
    447     for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
    448       int iColumn = column[j];
    449       index[numberNonZero++]=iColumn;
    450       array[iColumn] = value;
    451     }
    452     for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
    453       int iColumn = column[j];
    454       index[numberNonZero++]=iColumn;
    455       array[iColumn] = -value;
    456     }
    457     value = pi[iRow1]*scalar;
    458     for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
    459       int iColumn = column[j];
    460       double value2= array[iColumn];
    461       if (value2) {
    462         value2 += value;
    463       } else {
    464         value2 = value;
     676    int j;
     677    if (packed) {
     678      double pi0 = pi[0];
     679      double pi1 = pi[1];
     680      if (startPositive[iRow0+1]-startPositive[iRow0]>
     681          startPositive[iRow1+1]-startPositive[iRow1]) {
     682        int temp = iRow0;
     683        iRow0 = iRow1;
     684        iRow1 = temp;
     685        pi0=pi1;
     686        pi1=pi[0];
     687      }
     688      // and set up mark as char array
     689      char * marked = (char *) (index+columnArray->capacity());
     690      int * lookup = y->getIndices();
     691      double value = pi0*scalar;
     692      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
     693        int iColumn = column[j];
     694        array[numberNonZero] = value;
     695        marked[iColumn]=1;
     696        lookup[iColumn]=numberNonZero;
    465697        index[numberNonZero++]=iColumn;
    466698      }
    467       array[iColumn] = value2;
    468     }
    469     for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
    470       int iColumn = column[j];
    471       double value2= array[iColumn];
    472       if (value2) {
    473         value2 -= value;
    474       } else {
    475         value2 = -value;
     699      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
     700        int iColumn = column[j];
     701        array[numberNonZero] = -value;
     702        marked[iColumn]=1;
     703        lookup[iColumn]=numberNonZero;
    476704        index[numberNonZero++]=iColumn;
    477705      }
    478       array[iColumn] = value2;
    479     }
    480     // get rid of tiny values and zero out marked
    481     numberOriginal=numberNonZero;
    482     numberNonZero=0;
    483     for (i=0;i<numberOriginal;i++) {
    484       int iColumn = index[i];
    485       if (fabs(array[iColumn])>zeroTolerance) {
     706      int numberOriginal = numberNonZero;
     707      value = pi1*scalar;
     708      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
     709        int iColumn = column[j];
     710        if (marked[iColumn]) {
     711          int iLookup = lookup[iColumn];
     712          array[iLookup] += value;
     713        } else {
     714          if (fabs(value)>zeroTolerance) {
     715            array[numberNonZero] = value;
     716            index[numberNonZero++]=iColumn;
     717          }
     718        }
     719      }
     720      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
     721        int iColumn = column[j];
     722        if (marked[iColumn]) {
     723          int iLookup = lookup[iColumn];
     724          array[iLookup] -= value;
     725        } else {
     726          if (fabs(value)>zeroTolerance) {
     727            array[numberNonZero] = -value;
     728            index[numberNonZero++]=iColumn;
     729          }
     730        }
     731      }
     732      // get rid of tiny values and zero out marked
     733      int nDelete=0;
     734      for (j=0;j<numberOriginal;j++) {
     735        int iColumn = index[j];
     736        marked[iColumn]=0;
     737        if (fabs(array[j])<=zeroTolerance)
     738          nDelete++;
     739      }
     740      if (nDelete) {
     741        numberOriginal=numberNonZero;
     742        numberNonZero=0;
     743        for (j=0;j<numberOriginal;j++) {
     744          int iColumn = index[j];
     745          double value = array[j];
     746          array[j]=0.0;
     747          if (fabs(value)>zeroTolerance) {
     748            array[numberNonZero]=value;
     749            index[numberNonZero++]=iColumn;
     750          }
     751        }
     752      }
     753    } else {
     754      if (startPositive[iRow0+1]-startPositive[iRow0]<
     755          startPositive[iRow1+1]-startPositive[iRow1]) {
     756        int temp = iRow0;
     757        iRow0 = iRow1;
     758        iRow1 = temp;
     759      }
     760      int numberOriginal;
     761      int i;
     762      numberNonZero=0;
     763      double value;
     764      value = pi[iRow0]*scalar;
     765      CoinBigIndex j;
     766      for (j=startPositive[iRow0];j<startNegative[iRow0];j++) {
     767        int iColumn = column[j];
    486768        index[numberNonZero++]=iColumn;
    487       } else {
    488         array[iColumn]=0.0;
     769        array[iColumn] = value;
     770      }
     771      for (j=startNegative[iRow0];j<startPositive[iRow0+1];j++) {
     772        int iColumn = column[j];
     773        index[numberNonZero++]=iColumn;
     774        array[iColumn] = -value;
     775      }
     776      value = pi[iRow1]*scalar;
     777      for (j=startPositive[iRow1];j<startNegative[iRow1];j++) {
     778        int iColumn = column[j];
     779        double value2= array[iColumn];
     780        if (value2) {
     781          value2 += value;
     782        } else {
     783          value2 = value;
     784          index[numberNonZero++]=iColumn;
     785        }
     786        array[iColumn] = value2;
     787      }
     788      for (j=startNegative[iRow1];j<startPositive[iRow1+1];j++) {
     789        int iColumn = column[j];
     790        double value2= array[iColumn];
     791        if (value2) {
     792          value2 -= value;
     793        } else {
     794          value2 = -value;
     795          index[numberNonZero++]=iColumn;
     796        }
     797        array[iColumn] = value2;
     798      }
     799      // get rid of tiny values and zero out marked
     800      numberOriginal=numberNonZero;
     801      numberNonZero=0;
     802      for (i=0;i<numberOriginal;i++) {
     803        int iColumn = index[i];
     804        if (fabs(array[iColumn])>zeroTolerance) {
     805          index[numberNonZero++]=iColumn;
     806        } else {
     807          array[iColumn]=0.0;
     808        }
    489809      }
    490810    }
     
    495815    double value;
    496816    iRow = whichRow[0];
    497     value = pi[iRow]*scalar;
    498     if (fabs(value)>zeroTolerance) {
    499       CoinBigIndex j;
    500       for (j=startPositive[iRow];j<startNegative[iRow];j++) {
    501         int iColumn = column[j];
    502         index[numberNonZero++]=iColumn;
    503         array[iColumn] = value;
    504       }
    505       for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
    506         int iColumn = column[j];
    507         index[numberNonZero++]=iColumn;
    508         array[iColumn] = -value;
     817    CoinBigIndex j;
     818    if (packed) {
     819      value = pi[0]*scalar;
     820      if (fabs(value)>zeroTolerance) {
     821        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     822          int iColumn = column[j];
     823          array[numberNonZero] = value;
     824          index[numberNonZero++]=iColumn;
     825        }
     826        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     827          int iColumn = column[j];
     828          array[numberNonZero] = -value;
     829          index[numberNonZero++]=iColumn;
     830        }
     831      }
     832    } else {
     833      value = pi[iRow]*scalar;
     834      if (fabs(value)>zeroTolerance) {
     835        for (j=startPositive[iRow];j<startNegative[iRow];j++) {
     836          int iColumn = column[j];
     837          array[iColumn] = value;
     838          index[numberNonZero++]=iColumn;
     839        }
     840        for (j=startNegative[iRow];j<startPositive[iRow+1];j++) {
     841          int iColumn = column[j];
     842          array[iColumn] = -value;
     843          index[numberNonZero++]=iColumn;
     844        }
    509845      }
    510846    }
     
    532868  int numberToDo = y->getNumElements();
    533869  const int * which = y->getIndices();
    534   for (jColumn=0;jColumn<numberToDo;jColumn++) {
    535     int iColumn = which[jColumn];
    536     double value = 0.0;
    537     CoinBigIndex j=startPositive_[iColumn];
    538     for (;j<startNegative_[iColumn];j++) {
    539       int iRow = indices_[j];
    540       value += pi[iRow];
    541     }
    542     for (;j<startPositive_[iColumn+1];j++) {
    543       int iRow = indices_[j];
    544       value -= pi[iRow];
    545     }
    546     if (fabs(value)>zeroTolerance) {
    547       index[numberNonZero++]=iColumn;
    548       array[iColumn]=value;
     870  bool packed = rowArray->packedMode();
     871  if (packed) {
     872    // need to expand pi into y
     873    int numberInRowArray = rowArray->getNumElements();
     874    assert(y->capacity()>=model->numberRows());
     875    double * piOld = pi;
     876    pi = y->denseVector();
     877    const int * whichRow = rowArray->getIndices();
     878    int i;
     879    for (i=0;i<numberInRowArray;i++) {
     880      int iRow = whichRow[i];
     881      pi[iRow]=piOld[i];
     882    }
     883    // Must line up with y
     884    for (jColumn=0;jColumn<numberToDo;jColumn++) {
     885      int iColumn = which[jColumn];
     886      double value = 0.0;
     887      CoinBigIndex j=startPositive_[iColumn];
     888      for (;j<startNegative_[iColumn];j++) {
     889        int iRow = indices_[j];
     890        value += pi[iRow];
     891      }
     892      for (;j<startPositive_[iColumn+1];j++) {
     893        int iRow = indices_[j];
     894        value -= pi[iRow];
     895      }
     896      array[jColumn]=value;
     897    }
     898    for (i=0;i<numberInRowArray;i++) {
     899      int iRow = whichRow[i];
     900      pi[iRow]=0.0;
     901    }
     902  } else {
     903    for (jColumn=0;jColumn<numberToDo;jColumn++) {
     904      int iColumn = which[jColumn];
     905      double value = 0.0;
     906      CoinBigIndex j=startPositive_[iColumn];
     907      for (;j<startNegative_[iColumn];j++) {
     908        int iRow = indices_[j];
     909        value += pi[iRow];
     910      }
     911      for (;j<startPositive_[iColumn+1];j++) {
     912        int iRow = indices_[j];
     913        value -= pi[iRow];
     914      }
     915      if (fabs(value)>zeroTolerance) {
     916        index[numberNonZero++]=iColumn;
     917        array[iColumn]=value;
     918      }
    549919    }
    550920  }
     
    598968  return numberElements;
    599969}
     970/* If element NULL returns number of elements in column part of basis,
     971   If not NULL fills in as well */
     972CoinBigIndex
     973ClpPlusMinusOneMatrix::fillBasis(const ClpSimplex * model,
     974                                 const int * whichColumn,
     975                                 int numberBasic,
     976                                 int numberColumnBasic,
     977                                 int * indexRowU, int * indexColumnU,
     978                                 double * elementU) const
     979{
     980  int i;
     981  CoinBigIndex numberElements=0;
     982  if (elementU!=NULL) {
     983    assert (columnOrdered_);
     984    for (i=0;i<numberColumnBasic;i++) {
     985      int iColumn = whichColumn[i];
     986      CoinBigIndex j=startPositive_[iColumn];
     987      for (;j<startNegative_[iColumn];j++) {
     988        int iRow = indices_[j];
     989        indexRowU[numberElements]=iRow;
     990        indexColumnU[numberElements]=numberBasic;
     991        elementU[numberElements++]=1.0;
     992      }
     993      for (;j<startPositive_[iColumn+1];j++) {
     994        int iRow = indices_[j];
     995        indexRowU[numberElements]=iRow;
     996        indexColumnU[numberElements]=numberBasic;
     997        elementU[numberElements++]=-1.0;
     998      }
     999      numberBasic++;
     1000    }
     1001  } else {
     1002    for (i=0;i<numberColumnBasic;i++) {
     1003      int iColumn = whichColumn[i];
     1004      numberElements += startPositive_[iColumn+1]-startPositive_[iColumn];
     1005    }
     1006  }
     1007  return numberElements;
     1008}
    6001009/* Unpacks a column into an CoinIndexedvector
    601       Note that model is NOT const.  Bounds and objective could
    602       be modified if doing column generation */
     1010 */
    6031011void
    6041012ClpPlusMinusOneMatrix::unpack(const ClpSimplex * model,
     
    6151023    rowArray->add(iRow,-1.0);
    6161024  }
     1025}
     1026/* Unpacks a column into an CoinIndexedvector
     1027** in packed foramt
     1028Note that model is NOT const.  Bounds and objective could
     1029be modified if doing column generation (just for this variable) */
     1030void
     1031ClpPlusMinusOneMatrix::unpackPacked(ClpSimplex * model,
     1032                            CoinIndexedVector * rowArray,
     1033                            int iColumn) const
     1034{
     1035  int * index = rowArray->getIndices();
     1036  double * array = rowArray->denseVector();
     1037  int number = 0;
     1038  CoinBigIndex j=startPositive_[iColumn];
     1039  for (;j<startNegative_[iColumn];j++) {
     1040    int iRow = indices_[j];
     1041    array[number]=1.0;
     1042    index[number++]=iRow;
     1043  }
     1044  for (;j<startPositive_[iColumn+1];j++) {
     1045    int iRow = indices_[j];
     1046    array[number]=-1.0;
     1047    index[number++]=iRow;
     1048  }
     1049  rowArray->setNumElements(number);
     1050  rowArray->setPackedMode(true);
    6171051}
    6181052/* Adds multiple of a column into an CoinIndexedvector
     
    6931127ClpPlusMinusOneMatrix::deleteCols(const int numDel, const int * indDel)
    6941128{
    695   abort();
     1129  int iColumn;
     1130  int newSize=startPositive_[numberColumns_];;
     1131  int numberBad=0;
     1132  // Use array to make sure we can have duplicates
     1133  int * which = new int[numberColumns_];
     1134  memset(which,0,numberColumns_*sizeof(int));
     1135  int nDuplicate=0;
     1136  for (iColumn=0;iColumn<numDel;iColumn++) {
     1137    int jColumn = indDel[iColumn];
     1138    if (jColumn<0||jColumn>=numberColumns_) {
     1139      numberBad++;
     1140    } else {
     1141      newSize -= startPositive_[jColumn+1]-startPositive_[jColumn];
     1142      if (which[jColumn])
     1143        nDuplicate++;
     1144      else
     1145        which[jColumn]=1;
     1146    }
     1147  }
     1148  if (numberBad)
     1149    throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
     1150  int newNumber = numberColumns_-numDel+nDuplicate;
     1151  // Get rid of temporary arrays
     1152  delete [] lengths_;
     1153  lengths_=NULL;
     1154  delete [] elements_;
     1155  elements_= NULL;
     1156  int * newPositive = new int [newNumber+1];
     1157  int * newNegative = new int [newNumber];
     1158  int * newIndices = new int [newSize];
     1159  newNumber=0;
     1160  newSize=0;
     1161  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1162    if (!which[iColumn]) {
     1163      int start,end;
     1164      int i;
     1165      start = startPositive_[iColumn];
     1166      end=startNegative_[iColumn];
     1167      newPositive[newNumber]=newSize;
     1168      for (i=start;i<end;i++)
     1169        newIndices[newSize++]=indices_[i];
     1170      start = startNegative_[iColumn];
     1171      end=startPositive_[iColumn+1];
     1172      newNegative[newNumber++]=newSize;
     1173      for (i=start;i<end;i++)
     1174        newIndices[newSize++]=indices_[i];
     1175    }
     1176  }
     1177  newPositive[newNumber]=newSize;
     1178  delete [] which;
     1179  delete [] startPositive_;
     1180  startPositive_= newPositive;
     1181  delete [] startNegative_;
     1182  startNegative_= newNegative;
     1183  delete [] indices_;
     1184  indices_= newIndices;
     1185  numberColumns_ = newNumber;
    6961186}
    6971187/* Delete the rows whose indices are listed in <code>indDel</code>. */
     
    6991189ClpPlusMinusOneMatrix::deleteRows(const int numDel, const int * indDel)
    7001190{
    701   abort();
     1191  int iRow;
     1192  int numberBad=0;
     1193  // Use array to make sure we can have duplicates
     1194  int * which = new int[numberRows_];
     1195  memset(which,0,numberRows_*sizeof(int));
     1196  int nDuplicate=0;
     1197  for (iRow=0;iRow<numDel;iRow++) {
     1198    int jRow = indDel[iRow];
     1199    if (jRow<0||jRow>=numberRows_) {
     1200      numberBad++;
     1201    } else {
     1202      if (which[jRow])
     1203        nDuplicate++;
     1204      else
     1205        which[jRow]=1;
     1206    }
     1207  }
     1208  if (numberBad)
     1209    throw CoinError("Indices out of range", "deleteCols", "ClpPlusMinusOneMatrix");
     1210  int iElement;
     1211  int numberElements=startPositive_[numberColumns_];
     1212  int newSize=0;
     1213  for (iElement=0;iElement<numberElements;iElement++) {
     1214    iRow = indices_[iElement];
     1215    if (!which[iRow])
     1216      newSize++;
     1217  }
     1218  int newNumber = numberRows_-numDel+nDuplicate;
     1219  // Get rid of temporary arrays
     1220  delete [] lengths_;
     1221  lengths_=NULL;
     1222  delete [] elements_;
     1223  elements_= NULL;
     1224  int * newIndices = new int [newSize];
     1225  newSize=0;
     1226  int iColumn;
     1227  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1228    int start,end;
     1229    int i;
     1230    start = startPositive_[iColumn];
     1231    end=startNegative_[iColumn];
     1232    startPositive_[newNumber]=newSize;
     1233    for (i=start;i<end;i++) {
     1234        iRow = indices_[i];
     1235        if (!which[iRow])
     1236          newIndices[newSize++]=iRow;
     1237    }
     1238    start = startNegative_[iColumn];
     1239    end=startPositive_[iColumn+1];
     1240    startNegative_[newNumber]=newSize;
     1241    for (i=start;i<end;i++) {
     1242        iRow = indices_[i];
     1243        if (!which[iRow])
     1244          newIndices[newSize++]=iRow;
     1245    }
     1246  }
     1247  startPositive_[numberColumns_]=newSize;
     1248  delete [] which;
     1249  delete [] indices_;
     1250  indices_= newIndices;
     1251  numberRows_ = newNumber;
    7021252}
    7031253bool
     
    7291279  numberColumns_=numberColumns;
    7301280}
     1281/* Given positive integer weights for each row fills in sum of weights
     1282   for each column (and slack).
     1283   Returns weights vector
     1284*/
     1285CoinBigIndex *
     1286ClpPlusMinusOneMatrix::dubiousWeights(const ClpSimplex * model,int * inputWeights) const
     1287{
     1288  int numberRows = model->numberRows();
     1289  int numberColumns =model->numberColumns();
     1290  int number = numberRows+numberColumns;
     1291  CoinBigIndex * weights = new CoinBigIndex[number];
     1292  int i;
     1293  for (i=0;i<numberColumns;i++) {
     1294    CoinBigIndex j;
     1295    CoinBigIndex count=0;
     1296    for (j=startPositive_[i];j<startPositive_[i+1];j++) {
     1297      int iRow=indices_[j];
     1298      count += inputWeights[iRow];
     1299    }
     1300    weights[i]=count;
     1301  }
     1302  for (i=0;i<numberRows;i++) {
     1303    weights[i+numberColumns]=inputWeights[i];
     1304  }
     1305  return weights;
     1306}
     1307// Append Columns
     1308void
     1309ClpPlusMinusOneMatrix::appendCols(int number, const CoinPackedVectorBase * const * columns)
     1310{
     1311  int iColumn;
     1312  int size=0;
     1313  int numberBad=0;
     1314  for (iColumn=0;iColumn<number;iColumn++) {
     1315    int n=columns[iColumn]->getNumElements();
     1316    const double * element = columns[iColumn]->getElements();
     1317    size += n;
     1318    int i;
     1319    for (i=0;i<n;i++) {
     1320      if (fabs(element[i])!=1.0)
     1321        numberBad++;
     1322    }
     1323  }
     1324  if (numberBad)
     1325    throw CoinError("Not +- 1", "appendCols", "ClpPlusMinusOneMatrix");
     1326  // Get rid of temporary arrays
     1327  delete [] lengths_;
     1328  lengths_=NULL;
     1329  delete [] elements_;
     1330  elements_= NULL;
     1331  int numberNow = startPositive_[numberColumns_];
     1332  int * temp;
     1333  temp = new int [numberColumns_+1+number];
     1334  memcpy(temp,startPositive_,(numberColumns_+1)*sizeof(int));
     1335  delete [] startPositive_;
     1336  startPositive_= temp;
     1337  temp = new int [numberColumns_+number];
     1338  memcpy(temp,startNegative_,numberColumns_*sizeof(int));
     1339  delete [] startNegative_;
     1340  startNegative_= temp;
     1341  temp = new int [numberNow+size];
     1342  memcpy(temp,indices_,numberNow*sizeof(int));
     1343  delete [] indices_;
     1344  indices_= temp;
     1345  // now add
     1346  size=numberNow;
     1347  for (iColumn=0;iColumn<number;iColumn++) {
     1348    int n=columns[iColumn]->getNumElements();
     1349    const int * row = columns[iColumn]->getIndices();
     1350    const double * element = columns[iColumn]->getElements();
     1351    int i;
     1352    for (i=0;i<n;i++) {
     1353      if (element[i]==1.0)
     1354        indices_[size++]=row[i];
     1355    }
     1356    startNegative_[iColumn+numberColumns_]=size;
     1357    for (i=0;i<n;i++) {
     1358      if (element[i]==-1.0)
     1359        indices_[size++]=row[i];
     1360    }
     1361    startPositive_[iColumn+numberColumns_+1]=size;
     1362  }
     1363 
     1364  numberColumns_ += number;
     1365}
     1366// Append Rows
     1367void
     1368ClpPlusMinusOneMatrix::appendRows(int number, const CoinPackedVectorBase * const * rows)
     1369{
     1370  // Allocate arrays to use for counting
     1371  int * countPositive = new int [numberColumns_+1];
     1372  memset(countPositive,0,numberColumns_*sizeof(int));
     1373  int * countNegative = new int [numberColumns_];
     1374  memset(countNegative,0,numberColumns_*sizeof(int));
     1375  int iRow;
     1376  int size=0;
     1377  int numberBad=0;
     1378  for (iRow=0;iRow<number;iRow++) {
     1379    int n=rows[iRow]->getNumElements();
     1380    const int * column = rows[iRow]->getIndices();
     1381    const double * element = rows[iRow]->getElements();
     1382    size += n;
     1383    int i;
     1384    for (i=0;i<n;i++) {
     1385      int iColumn = column[i];
     1386      if (element[i]==1.0)
     1387        countPositive[iColumn++];
     1388      else if (element[i]==-1.0)
     1389        countNegative[iColumn++];
     1390      else
     1391        numberBad++;
     1392    }
     1393  }
     1394  if (numberBad)
     1395    throw CoinError("Not +- 1", "appendRows", "ClpPlusMinusOneMatrix");
     1396  // Get rid of temporary arrays
     1397  delete [] lengths_;
     1398  lengths_=NULL;
     1399  delete [] elements_;
     1400  elements_= NULL;
     1401  int numberNow = startPositive_[numberColumns_];
     1402  int * newIndices = new int [numberNow+size];
     1403  // Update starts and turn counts into positions
     1404  // also move current indices
     1405  int iColumn;
     1406  int numberAdded=0;
     1407  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1408    int n,now,move;
     1409    now = startPositive_[iColumn];
     1410    move = startNegative_[iColumn]-now;
     1411    n = countPositive[iColumn];
     1412    startPositive_[iColumn] += numberAdded;
     1413    memcpy(indices_+now,newIndices+startPositive_[iColumn],move*sizeof(int));
     1414    countPositive[iColumn]= startNegative_[iColumn]+numberAdded;
     1415    numberAdded += n;
     1416    now = startNegative_[iColumn];
     1417    move = startPositive_[iColumn+1]-now;
     1418    n = countNegative[iColumn];
     1419    startNegative_[iColumn] += numberAdded;
     1420    memcpy(indices_+now,newIndices+startNegative_[iColumn],move*sizeof(int));
     1421    countNegative[iColumn]= startPositive_[iColumn+1]+numberAdded;
     1422    numberAdded += n;
     1423  }
     1424  delete [] indices_;
     1425  indices_ = newIndices;
     1426  startPositive_[numberColumns_] += numberAdded;
     1427  // Now put in
     1428  for (iRow=0;iRow<number;iRow++) {
     1429    int newRow = numberRows_+iRow;
     1430    int n=rows[iRow]->getNumElements();
     1431    const int * column = rows[iRow]->getIndices();
     1432    const double * element = rows[iRow]->getElements();
     1433    int i;
     1434    for (i=0;i<n;i++) {
     1435      int iColumn = column[i];
     1436      int put;
     1437      if (element[i]==1.0) {
     1438        put = countPositive[iColumn];
     1439        countPositive[iColumn] = put+1;
     1440      } else {
     1441        put = countNegative[iColumn];
     1442        countNegative[iColumn] = put+1;
     1443      }
     1444      indices_[put]=newRow;
     1445    }
     1446  }
     1447  delete [] countPositive;
     1448  delete [] countNegative;
     1449  numberRows_ += number;
     1450}
     1451/* Returns largest and smallest elements of both signs.
     1452   Largest refers to largest absolute value.
     1453*/
     1454void
     1455ClpPlusMinusOneMatrix::rangeOfElements(double & smallestNegative, double & largestNegative,
     1456                       double & smallestPositive, double & largestPositive)
     1457{
     1458  int iColumn;
     1459  bool plusOne=false;
     1460  bool minusOne=false;
     1461  for (iColumn=0;iColumn<numberColumns_;iColumn++) {
     1462    if (startNegative_[iColumn]>startPositive_[iColumn])
     1463      plusOne=true;
     1464    if (startPositive_[iColumn+1]>startNegative_[iColumn])
     1465      minusOne=true;
     1466  }
     1467  if (minusOne) {
     1468    smallestNegative=-1.0;
     1469    largestNegative=-1.0;
     1470  } else {
     1471    smallestNegative=0.0;
     1472    largestNegative=0.0;
     1473  }
     1474  if (plusOne) {
     1475    smallestPositive=1.0;
     1476    largestPositive=1.0;
     1477  } else {
     1478    smallestPositive=0.0;
     1479    largestPositive=0.0;
     1480  }
     1481}
Note: See TracChangeset for help on using the changeset viewer.