C++ CoinPackedMatrix类代码示例

void OsiIF::loadDummyRow(OsiSolverInterface* s2, const double* lbounds, const double* ubounds, const double* objectives)
{
	CoinPackedVector *coinrow = new CoinPackedVector();
	CoinPackedMatrix *matrix =  new CoinPackedMatrix(false,0,0);
	matrix->setDimensions(0, numCols_);
	ArrayBuffer<int> dummy(1,false);
	dummy.push(0);
	char *senses = new char[1];
	double *rhs = new double[1];
	double *ranges = new double[1];
	coinrow->insert(0, 1.);
	matrix->appendRow(*coinrow);
	senses[0] = 'E';
	rhs[0] = 1.;
	ranges[0] = 0.;

	lpSolverTime_.start();
	s2->loadProblem(*matrix, lbounds, ubounds, objectives, senses, rhs, ranges);
	lpSolverTime_.stop();
	_remRows(dummy);

	delete coinrow;
	delete matrix;
	freeChar(senses);
	freeDouble(rhs);
	freeDouble(ranges);

	return;
}

/** Construct from a CoinPackedMatrix*/
TMat::TMat(const CoinPackedMatrix &M, MatrixStorageType T):
  iRow_(NULL), jCol_(NULL), value_(NULL), nnz_(M.getNumElements()),
  capacity_(M.getNumElements()), columnOrdering_(), rowOrdering_(),
  nonEmptyRows_(), nonEmptyCols_(){
  create(M);
  make_upper_triangular(T);
}

/* Load in an problem by copying the arguments (the constraints on the
   rows are given by lower and upper bounds). If a pointer is NULL then the
   following values are the default:
   <ul>
   <li> <code>colub</code>: all columns have upper bound infinity
   <li> <code>collb</code>: all columns have lower bound 0 
   <li> <code>rowub</code>: all rows have upper bound infinity
   <li> <code>rowlb</code>: all rows have lower bound -infinity
   <li> <code>obj</code>: all variables have 0 objective coefficient
      </ul>
*/
void 
CoinSnapshot::loadProblem(const CoinPackedMatrix& matrix,
			  const double* collb, const double* colub,   
			  const double* obj,
			  const double* rowlb, const double* rowub,
			  bool makeRowCopy)
{
  // Keep scalars (apart from objective value etc)
  gutsOfDestructor(3+8);
  numRows_ = matrix.getNumRows();
  numCols_ = matrix.getNumCols();
  numElements_ = matrix.getNumElements();
  owned_.matrixByCol = 1;
  matrixByCol_ = new CoinPackedMatrix(matrix);
  if (makeRowCopy) {
    owned_.matrixByRow = 1;
    CoinPackedMatrix * matrixByRow = new CoinPackedMatrix(matrix);
    matrixByRow->reverseOrdering();
    matrixByRow_ = matrixByRow;
  }
  colLower_ = CoinCopyOfArray(collb,numCols_,0.0);
  colUpper_ = CoinCopyOfArray(colub,numCols_,infinity_);
  objCoefficients_ = CoinCopyOfArray(obj,numCols_,0.0);
  rowLower_ = CoinCopyOfArray(rowlb,numRows_,-infinity_);
  rowUpper_ = CoinCopyOfArray(rowub,numRows_,infinity_);
  // do rhs as well
  createRightHandSide();
}

lpHook::lpHook(const double* clb, const double* cub, const double* obj,
               const double* rhs, const char* sense,
               const CoinPackedMatrix& mat)
{
     const int colnum = mat.getNumCols();
     const int rownum = mat.getNumRows();

     colupper_ = new double[colnum];
     collower_ = new double[colnum];
     objcoeffs_ = new double[colnum];
     rhs_ = new double[rownum];
     sense_ = new char[rownum];

     std::copy(clb, clb + colnum, collower_);
     std::copy(cub, cub + colnum, colupper_);
     std::copy(obj, obj + colnum, objcoeffs_);
     std::copy(rhs, rhs + rownum, rhs_);
     std::copy(sense, sense + rownum, sense_);

     if (mat.isColOrdered()) {
          colMatrix_.copyOf(mat);
          rowMatrix_.reverseOrderedCopyOf(mat);
     } else {
          rowMatrix_.copyOf(mat);
          colMatrix_.reverseOrderedCopyOf(mat);
     }
}

void
OsiVolSolverInterface::loadProblem(const CoinPackedMatrix& matrix,
				   const double* collb, const double* colub,
				   const double* obj,
				   const char* rowsen, const double* rowrhs,   
				   const double* rowrng)
{
   gutsOfDestructor_();
   const int rownum = matrix.getNumRows();
   const int colnum = matrix.getNumCols();

   if (matrix.isColOrdered()) {
      colMatrix_ = matrix;
      colMatrixCurrent_ = true;
      rowMatrixCurrent_ = false;
      maxNumcols_ = colMatrix_.getMaxMajorDim();
      maxNumrows_ = static_cast<int>((1+colMatrix_.getExtraGap()) *
				     colMatrix_.getMinorDim());
   } else {
      rowMatrix_ = matrix;
      rowMatrixCurrent_ = true;
      colMatrixCurrent_ = false;
      maxNumcols_ = static_cast<int>((1+rowMatrix_.getExtraGap()) *
				     rowMatrix_.getMinorDim());
      maxNumrows_ = rowMatrix_.getMaxMajorDim();
   }

   initFromRhsSenseRange(rownum, rowsen, rowrhs, rowrng);
   initFromClbCubObj(colnum, collb, colub, obj);
}

// take columns of matrix and add each to arrays for matrix under construction
void addSubMatr (int *start, int *len, int *ind, double *el,
		 CoinPackedMatrix &A, 
		 CoinPackedVector &v, 
		 int &cur,
		 int &ncols,
		 int dispM, 
		 int dispVec, 
		 int finalrow) {

  const int 
    *aLe  = A.getVectorLengths (),
    *aIn  = A.getIndices  (),
    *vIn  = v.getIndices  (),
     aCol = A.getMajorDim ();

  int vNum = v.getNumElements ();

  const double 
    *aEl = A.getElements (),
    *vEl = v.getElements ();

  // add each column
  for (int i=0; i<aCol; i++, len++) {

    *start++ = cur;
    *len     = *aLe++;

    // matrix entries
    for (int j = 0; j < *len; j++) {
      *ind++ = dispM + *aIn++;
      *el++  = *aEl++;
    }

    cur += *len;

    // check if there is a corresponding rhs
    if (vNum && (*vIn == i)) {

      ++*len;
      cur++;
      *ind++ = dispVec;
      *el++  = *vEl;

      vIn++;
      vEl++;
      --vNum;
    }

    // normalization entry
    ++*len;
    cur++;
    *ind++ = finalrow;
    *el++  = 1.;

    ++ncols;
  }

  *start = cur;
}

StochGenMatrix* sTreeImpl::createC() 
{
  //is this node a dead-end for this process?
  if(commWrkrs==MPI_COMM_NULL)
    return new StochGenDummyMatrix(m_id);

  StochGenMatrix* C = NULL;
  if (m_id==0) {
    CoinPackedMatrix Crow; 
    Crow.reverseOrderedCopyOf( in.getFirstStageConstraints() );

    // number of nz in the rows corresponding to ineq constraints
    int nnzD=countNNZ( Crow, 
		       in.getFirstStageRowLB(), 
		       in.getFirstStageRowUB(), 
		       ineq_comp());
    C = new StochGenMatrix( m_id, MZ, N, 
			    m_mz, 0,   0,    // C does not exist for the root
			    m_mz, m_nx, nnzD, // D is 1st stage ineq matrix
			    commWrkrs );
    extractRows( Crow,
		 in.getFirstStageRowLB(), in.getFirstStageRowUB(), ineq_comp(),
		 C->Bmat->krowM(), C->Bmat->jcolM(), C->Bmat->M() );
    //printf("  -- 1st stage mz=%lu nx=%lu nnzD=%d\n", m_mz, m_nx, nnzD);
  } else {
    int scen=m_id-1;
    CoinPackedMatrix Crow, Drow; 
    Crow.reverseOrderedCopyOf( in.getLinkingConstraints(scen) );
    Drow.reverseOrderedCopyOf( in.getSecondStageConstraints(scen) );

    int nnzC=countNNZ( Crow, in.getSecondStageRowLB(scen), 
		       in.getSecondStageRowUB(scen), ineq_comp() );
    int nnzD=countNNZ( Drow, in.getSecondStageRowLB(scen), 
		       in.getSecondStageRowUB(scen), ineq_comp() );

    C = new StochGenMatrix( m_id, MZ, N, 
			    m_mz, parent->m_nx, nnzC, 
			    m_mz, m_nx,         nnzD,
			    commWrkrs );
    //printf("  -- 2nd stage mz=%lu nx=%lu   1st stage nx=%lu     nnzC=%d nnzD=%d\n", m_mz, m_nx, parent->m_nx, nnzC, nnzD);
    extractRows( Crow,
		 in.getSecondStageRowLB(scen), 
		 in.getSecondStageRowUB(scen), 
		 ineq_comp(),
		 C->Amat->krowM(), C->Amat->jcolM(), C->Amat->M() );
    extractRows( Drow,
		 in.getSecondStageRowLB(scen), 
		 in.getSecondStageRowUB(scen), 
		 ineq_comp(),
		 C->Bmat->krowM(), C->Bmat->jcolM(), C->Bmat->M() );
  }
  
  for(size_t it=0; it<children.size(); it++) {
    StochGenMatrix* child = ((sTreeImpl*)children[it])->createC();
    C->AddChild(child);
  }
  return C;
}

// Create row-wise copy from MatrixByCol
void 
CoinSnapshot::createMatrixByRow()
{
  if (owned_.matrixByRow)
    delete matrixByRow_;
  assert (matrixByCol_);
  owned_.matrixByRow = 1;
  CoinPackedMatrix * matrixByRow = new CoinPackedMatrix(*matrixByCol_);
  matrixByRow->reverseOrdering();
  matrixByRow_ = matrixByRow;
}

int main () {
  // generate conic constraints
  CoinPackedMatrix * A = 0;
  CoinPackedVector * b = 0;
  CoinPackedVector * d = 0;
  double h = 0;
  generate_conic_constraint(A, b, d, h);
  A->dumpMatrix();
  // generate linear constraints
  CoinPackedMatrix * matrix = 0;
  double * col_lb = 0;
  double * col_ub = 0;
  double * obj = 0;
  double * row_lb = 0;
  double * row_ub = 0;
  generate_linear_constraints(matrix, col_lb, col_ub, obj, row_lb, row_ub);
  // create solver and load linear constraints
  OsiConicSolverInterface * solver = new ColaModel();
  // just add cols to the model
  solver->addCol(0, 0, 0, 0.0, solver->getInfinity(), 1.0, std::string("x1"));
  solver->addCol(0, 0, 0, 0.0, solver->getInfinity(), 0.0, std::string("x2"));
  solver->addCol(0, 0, 0, 0.0, solver->getInfinity(), 1.0, std::string("x3"));
  solver->addCol(0, 0, 0, 0.0, solver->getInfinity(), 0.0, std::string("x4"));
  //solver->loadProblem(*matrix, col_lb, NULL, obj, row_lb, row_ub);
  // add conic constraints
  solver->addConicConstraint(A, b, d, h);
  // todo(aykut) implement writing problem to mps files
  //solver->writeMps("ex");
  // solve problem
  solver->initialSolve();
  // print status
  std::cout << "is optimal: " << solver->isProvenOptimal() << std::endl;
  // print solution
  ColaModel * cm = dynamic_cast<ColaModel*>(solver);
  cm->print_stats();
  double const * sol = solver->getColSolution();
  for (int i=0; i<4; ++i) {
    std::cout << "x" << i << " " << sol[i] << std::endl;
  }
  delete A;
  delete b;
  delete d;
  delete matrix;
  delete solver;
  return 0;
}

bool
OsiTestSolverInterface::test_zero_one_minusone_(const CoinPackedMatrix& m) const
{
  const int vecnum = m.getMajorDim();
  const double* elem = m.getElements();
  const int* start = m.getVectorStarts();
  const int* length = m.getVectorLengths();
  int i, j;
  for (i = 0; i < vecnum; ++i) {
    for (j = start[i] + length[i] - 1; j >= start[i]; --j) {
      const double val = elem[j];
      if (val != 1.0 && val != 0.0 && val != -1.0) {
	return false;
      }
    }
  }
  return true;
}

   /** Equivalence.
       Two matrices are equivalent if they are both by rows or both by columns,
       they have the same dimensions, and each vector is equivalent. 
       In this method the FloatEqual function operator can be specified. 
   */
   template <class FloatEqual> bool 
   isEquivalent(const CoinPackedMatrix& rhs, const FloatEqual& eq) const
   {
      // Both must be column order or both row ordered and must be of same size
      if ((isColOrdered() ^ rhs.isColOrdered()) ||
	  (getNumCols() != rhs.getNumCols()) ||
	  (getNumRows() != rhs.getNumRows()) ||
	  (getNumElements() != rhs.getNumElements()))
	 return false;
     
      for (int i=getMajorDim()-1; i >= 0; --i) {
        CoinShallowPackedVector pv = getVector(i);
        CoinShallowPackedVector rhsPv = rhs.getVector(i);
        if ( !pv.isEquivalent(rhsPv,eq) )
          return false;
      }
      return true;
   }

int main(int argc, const char *argv[])
{
     ClpSimplex  model;
     int status;
     // Keep names
     if (argc < 2) {
          status = model.readMps("hello.mps", true);
     } else {
          status = model.readMps(argv[1], true);
     }
     if (status)
          exit(10);

     int numberColumns = model.numberColumns();
     int numberRows = model.numberRows();

     if (numberColumns > 80 || numberRows > 80) {
          printf("model too large\n");
          exit(11);
     }
     printf("This prints x wherever a non-zero element exists in the matrix.\n\n\n");

     char x[81];

     int iRow;
     // get row copy
     CoinPackedMatrix rowCopy = *model.matrix();
     rowCopy.reverseOrdering();
     const int * column = rowCopy.getIndices();
     const int * rowLength = rowCopy.getVectorLengths();
     const CoinBigIndex * rowStart = rowCopy.getVectorStarts();

     x[numberColumns] = '\0';
     for (iRow = 0; iRow < numberRows; iRow++) {
          memset(x, ' ', numberColumns);
          for (int k = rowStart[iRow]; k < rowStart[iRow] + rowLength[iRow]; k++) {
               int iColumn = column[k];
               x[iColumn] = 'x';
          }
          printf("%s\n", x);
     }
     printf("\n\n");
     return 0;
}

/* Modifies djs to allow for quadratic.
   returns quadratic offset */
CoinWorkDouble
ClpInterior::quadraticDjs(CoinWorkDouble * djRegion, const CoinWorkDouble * solution,
                          CoinWorkDouble scaleFactor)
{
     CoinWorkDouble quadraticOffset = 0.0;
#ifndef NO_RTTI
     ClpQuadraticObjective * quadraticObj = (dynamic_cast< ClpQuadraticObjective*>(objective_));
#else
     ClpQuadraticObjective * quadraticObj = NULL;
     if (objective_->type() == 2)
          quadraticObj = (static_cast< ClpQuadraticObjective*>(objective_));
#endif
     if (quadraticObj) {
          CoinPackedMatrix * quadratic = quadraticObj->quadraticObjective();
          const int * columnQuadratic = quadratic->getIndices();
          const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
          const int * columnQuadraticLength = quadratic->getVectorLengths();
          double * quadraticElement = quadratic->getMutableElements();
          int numberColumns = quadratic->getNumCols();
          for (int iColumn = 0; iColumn < numberColumns; iColumn++) {
               CoinWorkDouble value = 0.0;
               for (CoinBigIndex j = columnQuadraticStart[iColumn];
                         j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
                    int jColumn = columnQuadratic[j];
                    CoinWorkDouble valueJ = solution[jColumn];
                    CoinWorkDouble elementValue = quadraticElement[j];
                    //value += valueI*valueJ*elementValue;
                    value += valueJ * elementValue;
                    quadraticOffset += solution[iColumn] * valueJ * elementValue;
               }
               djRegion[iColumn] += scaleFactor * value;
          }
     }
     return quadraticOffset;
}

 //TODO: ugh force both row/col format?
 void loadDataAndSolution(const CoinPackedMatrix & rowMatrix,
                          const CoinPackedMatrix & colMatrix,
                          const double           * collb, 
                          const double           * colub,   
                          const double           * obj,
                          const double           * rowlb, 
                          const double           * rowub,
                          const char             * colType,
                          const double           * primalSol,
                          const double             infinity){
    data_->setMatrixByRow(&rowMatrix);
    data_->setMatrixByCol(&colMatrix);      
    data_->setNrow(rowMatrix.getNumRows());
    data_->setNcol(rowMatrix.getNumCols());
    data_->setColLower(collb);
    data_->setColUpper(colub);
    data_->setRowLower(rowlb);
    data_->setRowUpper(rowub);
    data_->setObj(obj);
    data_->setColType(colType);
    data_->setPrimalSol(primalSol);
    data_->setInfinity(infinity);
    data_->initializeOtherData();
 }

/* Orders rows and saves pointer to model */
int
ClpCholeskyWssmpKKT::order(ClpInterior * model)
{
    int numberRowsModel = model->numberRows();
    int numberColumns = model->numberColumns();
    int numberTotal = numberColumns + numberRowsModel;
    numberRows_ = 2 * numberRowsModel + numberColumns;
    rowsDropped_ = new char [numberRows_];
    memset(rowsDropped_, 0, numberRows_);
    numberRowsDropped_ = 0;
    model_ = model;
    CoinPackedMatrix * quadratic = NULL;
    ClpQuadraticObjective * quadraticObj =
        (dynamic_cast< ClpQuadraticObjective*>(model_->objectiveAsObject()));
    if (quadraticObj)
        quadratic = quadraticObj->quadraticObjective();
    int numberElements = model_->clpMatrix()->getNumElements();
    numberElements = numberElements + 2 * numberRowsModel + numberTotal;
    if (quadratic)
        numberElements += quadratic->getNumElements();
    // Space for starts
    choleskyStart_ = new CoinBigIndex[numberRows_+1];
    const CoinBigIndex * columnStart = model_->clpMatrix()->getVectorStarts();
    const int * columnLength = model_->clpMatrix()->getVectorLengths();
    const int * row = model_->clpMatrix()->getIndices();
    //const double * element = model_->clpMatrix()->getElements();
    // Now we have size - create arrays and fill in
    try {
        choleskyRow_ = new int [numberElements];
    } catch (...) {
        // no memory
        delete [] choleskyStart_;
        choleskyStart_ = NULL;
        return -1;
    }
    try {
        sparseFactor_ = new double[numberElements];
    } catch (...) {
        // no memory
        delete [] choleskyRow_;
        choleskyRow_ = NULL;
        delete [] choleskyStart_;
        choleskyStart_ = NULL;
        return -1;
    }
    int iRow, iColumn;

    sizeFactor_ = 0;
    // matrix
    if (!quadratic) {
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            choleskyStart_[iColumn] = sizeFactor_;
            choleskyRow_[sizeFactor_++] = iColumn;
            CoinBigIndex start = columnStart[iColumn];
            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
            for (CoinBigIndex j = start; j < end; j++) {
                choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
            }
        }
    } else {
        // Quadratic
        const int * columnQuadratic = quadratic->getIndices();
        const CoinBigIndex * columnQuadraticStart = quadratic->getVectorStarts();
        const int * columnQuadraticLength = quadratic->getVectorLengths();
        //const double * quadraticElement = quadratic->getElements();
        for (iColumn = 0; iColumn < numberColumns; iColumn++) {
            choleskyStart_[iColumn] = sizeFactor_;
            choleskyRow_[sizeFactor_++] = iColumn;
            for (CoinBigIndex j = columnQuadraticStart[iColumn];
                    j < columnQuadraticStart[iColumn] + columnQuadraticLength[iColumn]; j++) {
                int jColumn = columnQuadratic[j];
                if (jColumn > iColumn)
                    choleskyRow_[sizeFactor_++] = jColumn;
            }
            CoinBigIndex start = columnStart[iColumn];
            CoinBigIndex end = columnStart[iColumn] + columnLength[iColumn];
            for (CoinBigIndex j = start; j < end; j++) {
                choleskyRow_[sizeFactor_++] = row[j] + numberTotal;
            }
        }
    }
    // slacks
    for (; iColumn < numberTotal; iColumn++) {
        choleskyStart_[iColumn] = sizeFactor_;
        choleskyRow_[sizeFactor_++] = iColumn;
        choleskyRow_[sizeFactor_++] = iColumn - numberColumns + numberTotal;
    }
    // Transpose - nonzero diagonal (may regularize)
    for (iRow = 0; iRow < numberRowsModel; iRow++) {
        choleskyStart_[iRow+numberTotal] = sizeFactor_;
        // diagonal
        choleskyRow_[sizeFactor_++] = iRow + numberTotal;
    }
    choleskyStart_[numberRows_] = sizeFactor_;
    permuteInverse_ = new int [numberRows_];
    permute_ = new int[numberRows_];
    integerParameters_[0] = 0;
    int i0 = 0;
    int i1 = 1;
//.........这里部分代码省略.........