C++ BlockVector类代码示例

    // only need to accumulate blocks within number of circuit inputs
    // returns true if At is modified, false otherwise
    bool accumVector(BlockVector<T>& At) {
        const auto limit = std::min(At.stopIndex(), m_vec.size());
        if (At.startIndex() >= limit) {
            return false;

        } else {
            for (std::size_t i = At.startIndex(); i < limit; ++i) {
                m_vec[i] = At[3 + i] * m_random_rA;
#ifdef USE_ASSERT
                assert(! m_vec[i].isZero());
#endif
                At[3 + i] = T::zero();
            }

            return true;
        }
    }

  void computeh(double time, BlockVector& q0, SiconosVector& y)
  {

    std::cout <<"my_NewtonEulerR:: computeh" << std:: endl;
    std::cout <<"q0.size() = " << q0.size() << std:: endl;
    double height = q0.getValue(0) - _sBallRadius - q0.getValue(7);
    // std::cout <<"my_NewtonEulerR:: computeh _jachq" << std:: endl;
    // _jachq->display();

    y.setValue(0, height);
    _Nc->setValue(0, 1);
    _Nc->setValue(1, 0);
    _Nc->setValue(2, 0);
    _Pc1->setValue(0, q0.getValue(0) - _sBallRadius);
    _Pc1->setValue(1, q0.getValue(1));
    _Pc1->setValue(2, q0.getValue(2));

    _Pc2->setValue(0,q0.getValue(7));
    _Pc2->setValue(1,q0.getValue(8));
    _Pc2->setValue(2,q0.getValue(9));
    //printf("my_NewtonEulerR N, Pc\n");
    _Nc->display();
    _Pc1->display();
    _Pc2->display();
    std::cout <<"my_NewtonEulerR:: computeh ends" << std:: endl;
  }

  void computeh(double time, BlockVector& q0, SiconosVector& y)
  {
    double height = fabs(q0.getValue(0)) - _sBallRadius;
    // std::cout <<"my_NewtonEulerR:: computeh _jachq" << std:: endl;
    // _jachq->display();
    y.setValue(0, height);
    _Nc->setValue(0, 1);
    _Nc->setValue(1, 0);
    _Nc->setValue(2, 0);
    _Pc1->setValue(0, height);
    _Pc1->setValue(1, q0.getValue(1));
    _Pc1->setValue(2, q0.getValue(2));

    //_Pc2->setValue(0,hpc);
    //_Pc2->setValue(1,data[q0]->getValue(1));
    //_Pc2->setValue(2,data[q0]->getValue(2));
    //printf("my_NewtonEulerR N, Pc\n");
    //_Nc->display();
    //_Pc1->display();
  }

void prod(const SiconosMatrix& A, const BlockVector& x, SiconosVector& y, bool init)
{

  assert(!(A.isPLUFactorized()) && "A is PLUFactorized in prod !!");


  if (init)
    y.zero();
  unsigned int startRow = 0;
  unsigned int startCol = 0;
  // In private_addprod, the sum of all blocks of x, x[i], is computed: y = Sum_i (subA x[i]), with subA a submatrix of A,
  // starting from position startRow in rows and startCol in columns.
  // private_prod takes also into account the fact that each block of x can also be a block.
  VectorOfVectors::const_iterator it;
  for (it = x.begin(); it != x.end(); ++it)
  {
    private_addprod(A, startRow, startCol, **it, y);
    startCol += (*it)->size();
  }
}

void FirstOrderLinearR::computeg(double time, SiconosVector& lambda, SiconosVector& z, BlockVector& r)
{
  if (_pluginJacglambda->fPtr)
  {
    if (!_B)
      _B.reset(new SimpleMatrix(r.size(),lambda.size()));
    computeB(time, z, *_B);
  }

  prod(*_B, lambda, r, false);

}

void KneeJointR::computeh(double time, BlockVector& q0, SiconosVector& y)
{
  DEBUG_BEGIN("KneeJointR::computeh(double time, BlockVector& q0, SiconosVector& y)\n");
  DEBUG_EXPR(q0.display());

  double X1 = q0.getValue(0);
  double Y1 = q0.getValue(1);
  double Z1 = q0.getValue(2);
  double q10 = q0.getValue(3);
  double q11 = q0.getValue(4);
  double q12 = q0.getValue(5);
  double q13 = q0.getValue(6);
  DEBUG_PRINTF("X1 = %12.8e,\t Y1 = %12.8e,\t Z1 = %12.8e,\n",X1,Y1,Z1);
  DEBUG_PRINTF("q10 = %12.8e,\t q11 = %12.8e,\t q12 = %12.8e,\t q13 = %12.8e,\n",q10,q11,q12,q13);
  double X2 = 0;
  double Y2 = 0;
  double Z2 = 0;
  double q20 = 1;
  double q21 = 0;
  double q22 = 0;
  double q23 = 0;
  if(q0.getNumberOfBlocks()>1)
  {
    // SP::SiconosVector x2 = _d2->q();
    // DEBUG_EXPR( _d2->q()->display(););
    X2 = q0.getValue(7);
    Y2 = q0.getValue(8);
    Z2 = q0.getValue(9);
    q20 = q0.getValue(10);
    q21 = q0.getValue(11);
    q22 = q0.getValue(12);
    q23 = q0.getValue(13);
  }
  y.setValue(0, Hx(X1, Y1, Z1, q10, q11, q12, q13, X2, Y2, Z2, q20, q21, q22, q23));
  y.setValue(1, Hy(X1, Y1, Z1, q10, q11, q12, q13, X2, Y2, Z2, q20, q21, q22, q23));
  y.setValue(2, Hz(X1, Y1, Z1, q10, q11, q12, q13, X2, Y2, Z2, q20, q21, q22, q23));
  DEBUG_EXPR(y.display());
  DEBUG_END("KneeJointR::computeh(double time, BlockVector& q0, SiconosVector& y)\n");
    
}

// Copy from BlockVector
SiconosVector::SiconosVector(const BlockVector & vIn) : std11::enable_shared_from_this<SiconosVector>()
{
  if (ask<IsDense>(**(vIn.begin()))) // dense
  {
    _dense = true;
    vect.Dense = new DenseVect(vIn.size());
  }
  else
  {
    _dense = false;
    vect.Sparse = new SparseVect(vIn.size());
  }

  VectorOfVectors::const_iterator it;
  unsigned int pos = 0;
  for (it = vIn.begin(); it != vIn.end(); ++it)
  {
    setBlock(pos, **it);
    pos += (*it)->size();
  }

}

//-----------------------------------------------------------------------------
void BlockMatrix::mult(const BlockVector& x, BlockVector& y,
                       bool transposed) const
{
  if (transposed)
  {
    dolfin_error("BlockMatrix.cpp",
                 "compute transpose matrix-vector product",
                 "Not implemented for block matrices");
  }

  // Create temporary vector
  dolfin_assert(matrices[0][0]);

  // Loop over block rows
  for(std::size_t row = 0; row < matrices.shape()[0]; row++)
  {
    // RHS sub-vector
    GenericVector& _y = *(y.get_block(row));

    const GenericMatrix& _matA = *matrices[row][0];

    // Resize y and zero
    if (_y.empty())
      _matA.init_vector(_y, 0);
    _y.zero();

    // Loop over block columns
    std::shared_ptr<GenericVector>
      z_tmp = _matA.factory().create_vector();
    for(std::size_t col = 0; col < matrices.shape()[1]; ++col)
    {
      const GenericVector& _x = *(x.get_block(col));
      dolfin_assert(matrices[row][col]);
      matrices[row][col]->mult(_x, *z_tmp);
      _y += *z_tmp;
    }
  }
}

void FirstOrderLinearR::computeh(double time, BlockVector& x, SiconosVector& lambda,
                                 SiconosVector& z, SiconosVector& y)
{

  y.zero();

  if (_pluginJachx->fPtr)
  {
    if (!_C)
      _C.reset(new SimpleMatrix(y.size(),x.size()));
    computeC(time, z, *_C);
  }
  if (_pluginJachlambda->fPtr)
  {
    if (!_D)
      _D.reset(new SimpleMatrix(y.size(),lambda.size()));
    computeD(time, z, *_D);
  }
  if (_pluginf->fPtr)
  {
    if (!_F)
      _F.reset(new SimpleMatrix(y.size(),z.size()));
    computeF(time, z, *_F);
  }
  if (_plugine->fPtr)
  {
    if (!_e)
      _e.reset(new SiconosVector(y.size()));
    computee(time, z, *_e);
  }

  if (_C)
    prod(*_C, x, y, false);

  if (_D)
    prod(*_D, lambda, y, false);

  if (_e)
    y += *_e;

  if (_F)
    prod(*_F, z, y, false);

}

    void accumQuery(const BlockVector<G1>& query,
                    const BlockVector<Fr>& scalar,
                    ProgressCallback* callback = nullptr)
    {
#ifdef USE_ASSERT
        assert(query.space() == scalar.space() &&
               query.block() == scalar.block());
#endif

        accumQuery(query.vec(),
                   scalar.vec(),
                   callback);
    }

// Quantise a subband in in-place transform order
// This version of quantise_subbands assumes multiple quantisers per subband.
// It may be used for either quantising slices or for quantising subbands with codeblocks
const Array2D quantise_subbands(const Array2D& coefficients, const BlockVector& qIndices) {
  const Index transformHeight = coefficients.shape()[0];
  const Index transformWidth = coefficients.shape()[1];
  // TO DO: Check numberOfSubbands=3n+1 ?
  const int numberOfSubbands = qIndices.size();
  const int waveletDepth = (numberOfSubbands-1)/3;
  Index stride, offset; // stride is subsampling factor, offset is subsampling phase
  Array2D result(coefficients.ranges());

  // Create a view of the coefficients, representing the LL subband, quantise it,
  // then assign the result a view of the results array. This puts the quantised
  // LL subband into the result array in in-place transform order.
  // ArrayIndices2D objects specify the subset of array elements within a view,
  // that is they specify the subsampling factor and subsampling phase.
  stride = pow(2, waveletDepth);
  const ArrayIndices2D LLindices = // LLindices specifies the samples in the LL subband
    indices[Range(0,transformHeight,stride)][Range(0,transformWidth,stride)];
  result[LLindices] =
    quantise_LLSubband(coefficients[LLindices], qIndices[0]);

  // Next quantise the other subbands
  // Note: Level numbers go from zero for the lowest ("DC") frequencies to depth for
  // the high frequencies. This corresponds to the convention in the VC-2 specification.
  // Subands go from zero ("DC") to numberOfSubbands-1 for HH at the highest level
  for (char level=1, band=1; level<=waveletDepth; ++level) {
    stride = pow(2, waveletDepth+1-level);
    offset = stride/2;
    // Create a view of coefficients corresponding to a subband, then quantise it
    //Quantise HL subband
    const ArrayIndices2D HLindices = // HLindices specifies the samples in the HL subband
      indices[Range(0,transformHeight,stride)][Range(offset,transformWidth,stride)];
    result[HLindices] = quantise_block(coefficients[HLindices], qIndices[band++]);
    //Quantise LH subband
    const ArrayIndices2D LHindices = // LHindices specifies the samples in the LH subband
      indices[Range(offset,transformHeight,stride)][Range(0,transformWidth,stride)];
    result[LHindices] = quantise_block(coefficients[LHindices], qIndices[band++]);
    //Quantise HH subband
    const ArrayIndices2D HHindices = // HHindices specifies the samples in the HH subband
      indices[Range(offset,transformHeight,stride)][Range(offset,transformWidth,stride)];
    result[HHindices] = quantise_block(coefficients[HHindices], qIndices[band++]);
  }

  return result;
}

 PPZK_QueryHK(const BlockVector<Fr>& qap_query)
     : PPZK_QueryHK{qap_query.space(), qap_query.block()[0]}
 {}

//------------------------------------------------------------------------------
bool postdominatesAll(const BasicBlock *block, const BlockVector &blocks,
                      const PostDominatorTree *pdt) {
  return std::all_of(
      blocks.begin(), blocks.end(),
      [block, pdt](BasicBlock *iter) { return pdt->dominates(block, iter); });
}

// Template specialization for BasicBlock.
template <> void dumpVector(const BlockVector &toDump) {
  errs() << "Size: " << toDump.size() << "\n";
  for (auto element : toDump)
    errs() << element->getName() << " -- ";
  errs() << "\n";
}

/** Extract the master block from <code>problem</code>.
 * The constructor also sets up the solver for the newly created master
 * block. The master block can only be extracted if all sub-blocks have
 * already been extracted.
 * @param problem The problem from which to extract the master.
 * @param blocks  The sub blocks that have already been extracted.
 */
BendersOpt::Block::Block(Problem const *problem, BlockVector const &blocks)
   : env(), number(-1), vars(0), rows(0), cplex(0), cb(0)
{
   IloNumVarArray problemVars = problem->getVariables();
   IloRangeArray problemRanges = problem->getRows();

   IloExpr masterObj(env);
   IloNumVarArray masterVars(env);
   IloRangeArray masterRows(env);

   // Find columns that do not intersect block variables and
   // copy them to the master block.
   IdxMap idxMap;
   RowSet rowSet;
   for (IloInt j = 0; j < problemVars.getSize(); ++j) {
      IloNumVar x = problemVars[j];
      if ( problem->getBlock(x) < 0 ) {
         // Column is not in a block. Copy it to the master.
         IloNumVar v(env, x.getLB(), x.getUB(), x.getType(), x.getName());
         varMap.insert(VarMap::value_type(v, x));
         masterObj += problem->getObjCoef(x) * v;

         idxMap[x] = masterVars.getSize();
         masterVars.add(v);
      }
      else {
         // Column is in a block. Collect all rows that intersect
         // this column.
         RowSet const &intersected = problem->getIntersectedRows(x);
         for (RowSet::const_iterator it = intersected.begin();
              it != intersected.end(); ++it)
            rowSet.insert(*it);
         idxMap[x] = -1;
      }
   }

   // Pick up the rows that we need to copy.
   // These are the rows that are only intersected by master variables,
   // that is, the rows that are not in any block's rowset.
   for (IloInt i = 0; i < problemRanges.getSize(); ++i) {
      IloRange r = problemRanges[i];
      if ( rowSet.find(r) == rowSet.end() ) {
         IloRange masterRow(env, r.getLB(), r.getUB(), r.getName());
         IloExpr lhs(env);
         for (IloExpr::LinearIterator it = r.getLinearIterator(); it.ok(); ++it)
         {
            lhs += it.getCoef() * masterVars[idxMap[it.getVar()]];
         }
         masterRow.setExpr(lhs);
         masterRows.add(masterRow);
      }
   }

   // Adjust variable indices in blocks so that reference to variables
   // in the original problem become references to variables in the master.
   for (BlockVector::const_iterator b = blocks.begin(); b != blocks.end(); ++b) {
      for (std::vector<FixData>::iterator it = (*b)->fixed.begin(); it != (*b)->fixed.end(); ++it)
         it->col = idxMap[problemVars[it->col]];
   }

   // Create the eta variables, one for each block.
   // See the comments at the top of this file for details about the
   // eta variables.
   IloInt const firsteta = masterVars.getSize();
   for (BlockVector::size_type i = 0; i < blocks.size(); ++i) {
      std::stringstream s;
      s << "_eta" << i;
      IloNumVar eta(env, 0.0, IloInfinity, s.str().c_str());
      masterObj += eta;
      masterVars.add(eta);
   }

   // Create model and solver instance
   vars = masterVars;
   rows = masterRows;
   IloModel model(env);
   model.add(obj = IloObjective(env, masterObj, problem->getObjSense()));
   model.add(vars);
   model.add(rows);
   cplex = IloCplex(model);

   cplex.use(cb = new (env) LazyConstraintCallback(env, this, blocks,
                                              firsteta));

   for (IloExpr::LinearIterator it = obj.getLinearIterator(); it.ok(); ++it)
      objMap.insert(ObjMap::value_type(it.getVar(), it.getCoef()));
}