void PeriodicConditionLM2D2N::CalculateLocalSystem(MatrixType& rLeftHandSideMatrix, VectorType& rRightHandSideVector, ProcessInfo& rCurrentProcessInfo)
{
    if(rLeftHandSideMatrix.size1() != 6 || rLeftHandSideMatrix.size2() != 6) rLeftHandSideMatrix.resize(6, 6,false);
    noalias(rLeftHandSideMatrix) = ZeroMatrix(6, 6);

    if(rRightHandSideVector.size() != 6) rRightHandSideVector.resize(6, false);
	noalias( rRightHandSideVector ) = ZeroVector(6);

	// Nodal IDs = [slave ID, master ID]
	
	GeometryType& geom = GetGeometry();

	// get current values and form the system matrix

	Vector currentValues(6);
	currentValues(0) = geom[0].FastGetSolutionStepValue(DISPLACEMENT_X);
	currentValues(1) = geom[0].FastGetSolutionStepValue(DISPLACEMENT_Y);
	currentValues(2) = geom[1].FastGetSolutionStepValue(DISPLACEMENT_X);
	currentValues(3) = geom[1].FastGetSolutionStepValue(DISPLACEMENT_Y);
	currentValues(4) = geom[0].FastGetSolutionStepValue(DISPLACEMENT_LAGRANGE_X);
	currentValues(5) = geom[0].FastGetSolutionStepValue(DISPLACEMENT_LAGRANGE_Y);
	
	//KRATOS_WATCH(currentValues);

	rLeftHandSideMatrix(4,0) =  1.0;
	rLeftHandSideMatrix(4,2) = -1.0;
	rLeftHandSideMatrix(0,4) =  1.0;
	rLeftHandSideMatrix(2,4) = -1.0;
	rLeftHandSideMatrix(5,1) =  1.0;
	rLeftHandSideMatrix(5,3) = -1.0;
	rLeftHandSideMatrix(1,5) =  1.0;
	rLeftHandSideMatrix(3,5) = -1.0;

	// form residual

	noalias(rRightHandSideVector) -= prod( rLeftHandSideMatrix, currentValues );
}

void LinearSensor::initialize(const Model& m)
{
  // Call initialize of base class
  ControlSensor::initialize(m);

  // consistency checks
  if (!_matC)
  {
    RuntimeException::selfThrow("LinearSensor::initialize - no C matrix was given");
  }

  unsigned int colC = _matC->size(1);
  unsigned int rowC = _matC->size(0);
  // What happen here if we have more than one DS ?
  // This may be unlikely to happen.
  //  _DS = _model->nonSmoothDynamicalSystem()->dynamicalSystemNumber(0);
  if (colC != _DS->getN())
  {
    RuntimeException::selfThrow(" LinearSensor::initialize - The number of column of the C matrix must be equal to the length of x");
  }
  if (_matD)
  {
    unsigned int rowD = _matD->size(0);
    if (rowC != rowD)
    {
      RuntimeException::selfThrow("C and D must have the same number of rows");
    }
  }

  // --- Get the values ---
  // -> saved in a matrix data
  // -> event
  _storedY.reset(new SiconosVector(rowC));
  //  (_data[_eSensor])["StoredY"] = storedY;
  // set the dimension of the output
  *_storedY = prod(*_matC, *_DSx);
}

void CState::makeBoxes(const vec3 systemSize, const double interactionLength)
{
    ivec3 boxIndex;
    vec3 boxPos;

    NBoxes = conv_to<ivec>::from(floor(systemSize/interactionLength));
    boxDimensions = systemSize/NBoxes;
    nBoxes = prod(NBoxes);

    boxes = vector<CBox*>();
    boxes.reserve(nBoxes);
    boxes.resize(nBoxes); // using resize since we aren't using push_back to add items
    for (int ix = 0; ix < NBoxes(0); ix++)
    {
        for (int iy = 0; iy < NBoxes(1); iy++)
        {
            for (int iz = 0; iz < NBoxes(2); iz++)
            {
                boxIndex << ix << iy << iz;
                boxPos = boxIndex%boxDimensions;
                CBox* box = new CBox(nBoxes, boxPos, boxDimensions, boxIndex, NBoxes);
                boxes[calculate_box_number(boxIndex, NBoxes)] = box;
            }
        }
    }

    // finding all the neighbouring boxes of each box, taking into account the
    // periodic boundary conditions
    for (int i = 0; i < nBoxes; i++)
        boxes[i]->findNeighbours(systemSize);

    cout << "box dimensions (MD): " << boxDimensions.t()
         << "               (SI): " << boxDimensions.t()*L0
         << "system size    (MD): " << systemSize.t()
         << "               (SI): " << systemSize.t()*L0
         << "number of boxes    :    " << nBoxes << endl << endl;
}

void CombatCamera::ViewportRay(double screen_x, double screen_y,
                               double out_ray_origin[3], double out_ray_direction[3]) const
{
    Matrix projection(4, 4);
    std::copy(m_camera.getProjectionMatrix()[0], m_camera.getProjectionMatrix()[0] + 16,
              projection.data().begin());

    Matrix view(4, 4);
    std::copy(m_camera.getViewMatrix(true)[0], m_camera.getViewMatrix(true)[0] + 16, view.data().begin());

    Matrix inverse_vp = Inverse4(prod(projection, view));

    double nx = (2.0 * screen_x) - 1.0;
    double ny = 1.0 - (2.0 * screen_y);
    Matrix near_point(3, 1);
    near_point(0, 0) = nx;
    near_point(1, 0) = ny;
    near_point(2, 0) = -1.0;
    // Use mid_point rather than far point to avoid issues with infinite projection
    Matrix mid_point(3, 1);
    mid_point(0, 0) = nx;
    mid_point(1, 0) = ny;
    mid_point(2, 0) = 0.0;

    // Get ray origin and ray target on near plane in world space
    Matrix ray_origin = Matrix4xVector3(inverse_vp, near_point);
    Matrix ray_target = Matrix4xVector3(inverse_vp, mid_point);

    Matrix ray_direction = ray_target - ray_origin;
    ray_direction /=
        ray_direction(0, 0) * ray_direction(0, 0) +
        ray_direction(1, 0) * ray_direction(1, 0) +
        ray_direction(2, 0) * ray_direction(2, 0);

    std::copy(ray_origin.data().begin(), ray_origin.data().end(), out_ray_origin);
    std::copy(ray_direction.data().begin(), ray_direction.data().end(), out_ray_direction);
}

// Update inverse a1 after row in matrix a has changed
// get new row out of array1 newRow
//
// a1 = old inverse
// newRow  = new row lRow in new matrix a_new
// returns Det(a_old)/Det(a_new)
doublevar InverseUpdateRow(Array2 <doublevar> & a1, const Array1 <doublevar> & newRow,
                           const int lRow, const int n)
{
  Array1 <doublevar> & tmpColL(tmp11);
  tmpColL.Resize(n);
  Array1 <doublevar> & prod(tmp12);
  prod.Resize(n);

  doublevar f=0.0;
  for(int i=0;i<n;++i) {
    f += newRow[i]*a1(i,lRow);
  }
  f =-1.0/f;

  for(int j=0;j<n;++j){
    prod[j]   =0.0;
    tmpColL[j]=a1(j,lRow);
    for(int i=0;i<n;++i) {
      prod[j] += newRow[i]*a1(i,j);
    }
    prod[j] *= f;
  }

  for(int i=0;i<n;++i) {
    doublevar & t(tmpColL[i]);
    for(int j=0;j<n;++j) {
      a1(i,j) += t*prod[j];
    }
  }

  f = -f;
  for(int i=0;i<n;++i) {
    a1(i,lRow) = f*tmpColL[i];
  }
  return f;
}

// This algorithm factors each element a ∈ S over P if P is a base for S; otherwise it proclaims failure.
//
// Algorithm 21.2  [PDF page 27](http://cr.yp.to/lineartime/dcba-20040404.pdf)
//
// See [findfactors test](test-findfactors.html) for basic usage.
void find_factors(mpz_array *out, mpz_t *s, size_t from, size_t to, mpz_array *p) {
	mpz_t x, y, z;
	mpz_array d, q;
	size_t i, n = to - from;

	mpz_init(x);
	array_prod(p, x);

	mpz_init(y);
	prod(y, s, from, to);

	mpz_init(z);
	ppi(z, x, y);

	array_init(&d, p->size);
	array_split(&d, z, p);

	array_init(&q, p->size);
	for (i = 0; i < p->used; i++) {
		if (mpz_cmp(d.array[i], p->array[i]) == 0)
			array_add(&q, p->array[i]);
	}

	if (n == 0) {
		array_find_factor(out, y, &q);
	} else {
		find_factors(out, s, from, to - n/2 - 1, &q);
		find_factors(out, s, to - n/2, to, &q);
	}

	mpz_clear(x);
	mpz_clear(y);
	mpz_clear(z);
	array_clear(&d);
	array_clear(&q);
}

  ProbabilityDistribution* 
  StudentTProcessJeffreys::prediction(const vectord &query )
  {
    clock_t start = clock();
    double kq = computeSelfCorrelation(query);
    //    vectord kn = computeCrossCorrelation(query);
    mKernel.computeCrossCorrelation(mData.mX,query,mKn);
    vectord phi = mMean.getFeatures(query);
  
    //    vectord v(mKn);
    inplace_solve(mL,mKn,ublas::lower_tag());

    vectord rho = phi - prod(mKn,mKF);

    //    vectord rho(rq);
    inplace_solve(mL2,rho,ublas::lower_tag());
    
    double yPred = inner_prod(phi,mWML) + inner_prod(mKn,mAlphaF);
    double sPred = sqrt( mSigma * (kq - inner_prod(mKn,mKn) 
				   + inner_prod(rho,rho)));

    d_->setMeanAndStd(yPred,sPred);
    return d_;
  }

void AbstractMaterialLaw<DIM>::ComputeCauchyStress(c_matrix<double,DIM,DIM>& rF,
                                                   double pressure,
                                                   c_matrix<double,DIM,DIM>& rSigma)
{
    double detF = Determinant(rF);

    c_matrix<double,DIM,DIM> C = prod(trans(rF), rF);
    c_matrix<double,DIM,DIM> invC = Inverse(C);

    c_matrix<double,DIM,DIM> T;

    static FourthOrderTensor<DIM,DIM,DIM,DIM> dTdE; // not filled in, made static for efficiency

    ComputeStressAndStressDerivative(C, invC, pressure, T, dTdE, false);

    /*
     * Looping is probably more eficient then doing rSigma = (1/detF)*rF*T*transpose(rF),
     * which doesn't seem to compile anyway, as rF is a Tensor<2,DIM> and T is a
     * SymmetricTensor<2,DIM>.
     */
    for (unsigned i=0; i<DIM; i++)
    {
        for (unsigned j=0; j<DIM; j++)
        {
            rSigma(i,j) = 0.0;
            for (unsigned M=0; M<DIM; M++)
            {
                for (unsigned N=0; N<DIM; N++)
                {
                    rSigma(i,j) += rF(i,M)*T(M,N)*rF(j,N);
                }
            }
            rSigma(i,j) /= detF;
        }
    }
}

void SlidingReducedOrderObserver::process()
{
  DEBUG_BEGIN("void SlidingReducedOrderObserver::process()\n");
  if (!_pass)
  {
    DEBUG_PRINT("First pass \n ");
    _pass = true;
    //update the estimate using the first value of y, such that C\hat{x}_0 = y_0
    const SiconosVector& y = _sensor->y();
    _e->zero();
    prod(*_C, *_xHat, *_e);
    *_e -= y;

    SiconosVector tmpV(_DS->n());
    SimpleMatrix tmpC(*_C);
    for (unsigned int i = 0; i < _e->size(); ++i)
      tmpV(i) = (*_e)(i);

    tmpC.SolveByLeastSquares(tmpV);
    *(_xHat) -= tmpV;
    *(_DS->x()) -= tmpV;
    _DS->initMemory(1);
    _DS->swapInMemory();
    DEBUG_EXPR(_DS->display(););

microseconds NestedExpr::boost_ublas(const Args & args, std::mt19937 & gen)
{
	std::cout << "Test: boost uBLAS ";
	const NestedExprArgs & cur_args = dynamic_cast<const NestedExprArgs&>(args);

	uint32_t m = cur_args.matrix_size;
	uint32_t matr_size = m * m;

	boost::numeric::ublas::matrix<double> A(m, m), B(m, m), C(m, m), D(m, m), E(m,m);

	initialize(A.data().begin(), A.data().end(), gen);
	initialize(B.data().begin(), B.data().end(), gen);
	initialize(C.data().begin(), C.data().end(), gen);
	initialize(D.data().begin(), D.data().end(), gen);

	auto start = std::chrono::high_resolution_clock::now();
	noalias(E) = prod( A + B, C - D );
	auto end = std::chrono::high_resolution_clock::now();

	auto time = std::chrono::duration_cast<std::chrono::microseconds>( end - start);
	std::cout << time.count() << std::endl;

	return time;
}

static void fftn_c2r(const Rcomplex *z, R_len_t rank, const R_len_t *N,
                     double *res) {
  SEXP rTrue, cA, dim, Res;
  R_len_t n = prod(rank, N), i;

  rTrue = PROTECT(allocVector(LGLSXP, 1));
  LOGICAL(rTrue)[0] = 1;

  cA = PROTECT(allocVector(CPLXSXP, n));
  memcpy(COMPLEX(cA), z, sizeof(Rcomplex) * n);

  dim = PROTECT(allocVector(INTSXP, rank));
  memcpy(INTEGER(dim), N, sizeof(R_len_t) * rank);
  setAttrib(cA, R_DimSymbol, dim);

  Res = PROTECT(eval(lang3(install("fft"), cA, rTrue), R_GlobalEnv));

  /* Return result */
  for (i = 0; i < n; ++i)
    res[i] = COMPLEX(Res)[i].r;

  /* Unprotect all */
  UNPROTECT(4);
}

int vp(const my_point &a, const my_point &b, const my_point &c)
{
    double t[12];

    std::pair<double, double> tmp;
    tmp = prod(b.x, c.y);
    t[0] = tmp.first;
    t[1] = tmp.second;
    tmp = prod(-b.x, a.y);
    t[2] = tmp.first;
    t[3] = tmp.second;
    tmp = prod(-a.x, c.y);
    t[4] = tmp.first;
    t[5] = tmp.second;
    tmp = prod(-b.y, c.x);
    t[6] = tmp.first;
    t[7] = tmp.second;
    tmp = prod(b.y, a.x);
    t[8] = tmp.first;
    t[9] = tmp.second;
    tmp = prod(a.y, c.x);
    t[10] = tmp.first;
    t[11] = tmp.second;

    exp_sum(t, t + 2, 2, 2);
    exp_sum(t + 4, t + 6, 2, 2);
    exp_sum(t + 8, t + 10, 2, 2);
    exp_sum(t, t + 4, 4, 4);
    exp_sum(t, t + 8, 8, 4);

    for (int i = 11; i >= 0; i--)
    {
        if (t[i] > 0)
            return 1;
        if (t[i] < 0)
            return -1;
    }
    return 0;
}

void PlaneStressJ2::CalculateMaterialResponse(const Vector& StrainVector,
        const Matrix& DeformationGradient,
        Vector& StressVector,
        Matrix& AlgorithmicTangent,
        const ProcessInfo& CurrentProcessInfo,
        const Properties& props,
        const GeometryType& geom,
        const Vector& ShapeFunctionsValues,
        bool CalculateStresses,
        int CalculateTangent,
        bool SaveInternalVariables)
{
    KRATOS_TRY

    mE = props.GetValue(YOUNG_MODULUS);
    mNU = props.GetValue(POISSON_RATIO);
    msigma_y = props.GetValue(YIELD_STRESS);
    mH = 0.0;
    mtheta = 0.0;

//        double theta = 0.0;

    //resize output quantities
    if (StressVector.size() != 3 && CalculateStresses == true) StressVector.resize(3, false);
    if (AlgorithmicTangent.size1() != 3 && CalculateTangent != 0) AlgorithmicTangent.resize(3, 3, false);


    array_1d<double, 3 > elastic_strain;
    noalias(elastic_strain) = StrainVector;
    noalias(elastic_strain) -= mOldPlasticStrain;
    //        KRATOS_WATCH(StrainVector);
    //        KRATOS_WATCH(mOldPlasticStrain);

    boost::numeric::ublas::bounded_matrix<double, 3, 3 > C, Cinv, P;
    CalculateElasticMatrix(C);
    CalculateInverseElasticMatrix(Cinv);
    CalculateP(P);

//                        KRATOS_WATCH(C);

    array_1d<double, 3 > s_trial = prod(C, elastic_strain);
    array_1d<double, 3 > xi_trial = s_trial;
    noalias(s_trial) -= mbeta_old;

    //                KRATOS_WATCH(compute_f_trial(xi_trial, malpha_old));

    //        double fbar2 = fbar_2(0.0, xi_trial);
    //        double fbar_value = sqrt(fbar_2(0.0, xi_trial));
    //        double r2_value = R_2(0.0, fbar_value, malpha_old);
    //        double aaa = fbar_value - sqrt(r2_value);
    //                KRATOS_WATCH(sqrt(r2_value));
    //                KRATOS_WATCH(fbar_value);
    //                KRATOS_WATCH(sqrt(r2_value));
    //                        KRATOS_WATCH(aaa);

    double H1 = (1.0 - mtheta) * mH;
    ////        KRATOS_WATCH(xi_trial);
    //                KRATOS_WATCH(mbeta_old)
    if (compute_f_trial(xi_trial, malpha_old) < 0) //elastic case
    {
        if (CalculateStresses == true) noalias(StressVector) = s_trial;
        if (CalculateTangent != 0) noalias(AlgorithmicTangent) = C;
        //note that in this case internal variables are not modified!

    }
    else
    {
        //algorithm copied identically from the Simo Hughes, BOX3.3, pag 130
        double dgamma = ComputeDGamma(xi_trial, malpha_old);
        double ccc = 0.6666666666666667 * dgamma*H1;

        //            KRATOS_WATCH(dgamma);
        //            KRATOS_WATCH(xi_trial);
        //                        KRATOS_WATCH(malpha_old);

        //calculate XImat
        //note that here i reuse the C as auxiliary variable as i don't need it anymore
        boost::numeric::ublas::bounded_matrix<double, 3, 3 > XImat;
        noalias(C) = Cinv;
        noalias(C) += (dgamma / (1.0 + ccc)) * P;
        double detC;
        InvertMatrix(C, XImat, detC);
        //            KRATOS_WATCH(XImat);


        array_1d<double, 3 > aux, xi;
        noalias(aux) = prod(Cinv, xi_trial);
        noalias(xi) = prod(XImat, aux);
        xi /= (1.0 + ccc);

        //            KRATOS_WATCH(compute_f_trial(xi, malpha_old));

        noalias(mbeta_n1) = mbeta_old;
        noalias(mbeta_n1) += ccc*xi;

        array_1d<double, 3 > stress;
        noalias(stress) = xi;
        noalias(stress) += mbeta_n1;
        if (CalculateStresses == true) noalias(StressVector) = s_trial;

//.........这里部分代码省略.........

int test_main (int, char *[])
{
    prod(matrix_3x1(), matrix_3x1());

    return 0;
}

//.........这里部分代码省略.........


	  //Initialize water capacity  matrix
	  ublas::banded_matrix<double> Cw (cell_size, cell_size, 0, 0);
	  for (size_t c = 0; c < cell_size; c++)
	    Cw (c, c) = soil.Cw2 (c, h[c]);
	  
          std::vector<double> h_std (cell_size);
          //ublas vector -> std vector 
          std::copy(h.begin (), h.end (), h_std.begin ());

#ifdef TEST_OM_DEN_ER_BRUGT
          for (size_t cell = 0; cell != cell_size ; ++cell) 
            {				
              S_macro (cell) = (S_matrix[cell] + S_drain[cell]) 
                * geo.cell_volume (cell);
            }
#endif

	  //Initialize sum matrix
	  Solver::Matrix summat (cell_size);  
	  noalias (summat) = diff + Dm_mat;

	  //Initialize sum vector
	  ublas::vector<double> sumvec (cell_size);  
	  sumvec = grav + B + Gm + Dm_vec - S_vol
#ifdef TEST_OM_DEN_ER_BRUGT
            - S_macro
#endif
            ; 

	  // QCw is shorthand for Qmatrix * Cw
	  Solver::Matrix Q_Cw (cell_size);
	  noalias (Q_Cw) = prod (Qmat, Cw);

	  //Initialize A-matrix
	  Solver::Matrix A (cell_size);  
	  noalias (A) = (1.0 / ddt) * Q_Cw - summat;  

	  // Q_Cw_h is shorthand for Qmatrix * Cw * h
	  const ublas::vector<double> Q_Cw_h = prod (Q_Cw, h);

	  //Initialize b-vector
	  ublas::vector<double> b (cell_size);  
	  //b = sumvec + (1.0 / ddt) * (Qmatrix * Cw * h + Qmatrix *(Wxx-Wyy));
	  b = sumvec + (1.0 / ddt) * (Q_Cw_h
				      + prod (Qmat, Theta_previous-Theta));

	  // Force active drains to zero h.
          drain (geo, drain_cell, drain_water_level,
		 h, Theta_previous, Theta, S_vol,
#ifdef TEST_OM_DEN_ER_BRUGT
                 S_macro,
#endif
                 dq, ddt, drain_cell_on, A, b, debug, msg);  
          
          try {
            solver->solve (A, b, h); // Solve Ah=b with regard to h.
          } catch (const char *const error) {
              std::ostringstream tmp;
              tmp << "Could not solve equation system: " << error;
              msg.warning (tmp.str ());
              // Try smaller timestep.
              iterations_used = max_loop_iter + 100;
              break;
          }

double ModelFunction::operator()(const vector< double > &x) const
{
  // NOTE: the constant is not included
  return inner_prod(x, prod(A_, x))/2.0 + inner_prod(b_, x);
}

const vector< double > &ModelFunction::g(const vector< double > &x, vector< double > &g)
{
  g = prod(A_, x) + b_;
  return g;
}

  Interval LDB::operator ()(const LDB::ipolynomial_type    *poly_ptr,
                            const LDB::additional_var_data *data_ptr  ) const
  {
    unsigned int data_size = data_ptr->size();

    /*
      Abspalten der linearen Terme. Die Konstante wird den hoeheren Termen zugerechnet.
      Gleichzeitig werden die hoeheren Terme eingeschlossen.
    */

    Interval rest(0.0);
    std::vector<Interval> arguments( data_ptr->size() );
    for(unsigned int i = 0; i < data_size; i++)
      if( (*data_ptr)[ i ].operator ->() != 0 )
        arguments[ i ] = ((*data_ptr)[ i ])->domain_ - ((*data_ptr)[ i ])->devel_point_;

    std::vector<Interval> linear_coeffs(data_size, Interval(0.0));
    std::vector<unsigned int> var_codes; var_codes.reserve(data_size); //Codes of Vars with linear term.

    LDB::ipolynomial_type::const_iterator
      curr = poly_ptr->begin(),
      last = poly_ptr->end();

    while( curr != last )
    {
      if( curr->key().expnt_sum() == 0 )      //Konstanter Term.
      {
        rest += curr->value();
      }
      else if( curr->key().expnt_sum() == 1 ) //Linearer Term.
      {
        unsigned int i = curr->key().min_var_code();

        var_codes.push_back(i-1);
        linear_coeffs[i-1] = curr->value();
      }
      else //Nonlinear term.
      {
        Interval prod(1.0);

        //Evaluate the current monomial.
        for(unsigned int i = curr->key().min_var_code(); i <= curr->key().max_var_code(); i++)
        {
          unsigned int e = curr->key().expnt_of_var(i);
          if( e != 0 )
          {
            prod *= power( arguments[i-1], e );
          }
        }
        //Multiply with coefficient.
        prod *= curr->value();

        //Add enclosure to bound interval.
        rest += prod;
      }

      ++curr;
    }

    if( var_codes.size() == 0 ) return rest; //No linear terms in polynomial. So return with enclosure of
    //higher order terms.

    /*
      Lokale Kopien der 'Domains' anlegen, da diese im weiteren Verlauf unter Umstaenden
      veraendert werden.
    */

    std::vector<Interval> domain_of_max(data_size,Interval(0.0));
    std::vector<Interval> domain_of_min(data_size,Interval(0.0));

    for(unsigned int i = 0; i < data_size; i++)
    {
      if( (*data_ptr)[i].operator ->() ) //There is information.
      {
        domain_of_max[i] = domain_of_min[i] = ((*data_ptr)[i])->domain_;
      }
    }

    /*
      Jetzt beginnt der Algorithmus des LDB range bounders.
    */

    Interval idelta(rest.diam());
    for(unsigned int i = 0; i < var_codes.size(); i++)
    {
      unsigned int k = var_codes[i];
      Interval b = linear_coeffs[k];
      idelta += b.diam() * abs( ((*data_ptr)[k])->domain_ - ((*data_ptr)[k])->devel_point_ );
    }
    double delta = idelta.sup();

    bool resizing = false;
    for(unsigned int i = 0; i < var_codes.size(); i++)
    {
      unsigned int k = var_codes[i];

      double   mid_b  = linear_coeffs[k].mid();
      Interval abs_b  = Interval( std::abs(mid_b) );
      Interval domain = ((*data_ptr)[k])->domain_;
      Interval upper  = abs_b * domain.diam();
//.........这里部分代码省略.........

 types::ndarray<T,N - 1>
 prod(types::ndarray<T,N> const& array, long axis)
 {
     if(axis<0 || axis >=long(N))
         throw types::ValueError("axis out of bounds");
     auto shape = array.shape;
     if(axis==0)
     {
         types::array<long, N> shp;
         shp[0] = 1;
         std::copy(shape.begin() + 1, shape.end(), shp.begin() + 1);
         types::ndarray<T,N> out(shp, 1);
         return std::accumulate(array.begin(), array.end(), *out.begin(), proxy::multiply());
     }
     else
     {
         types::array<long, N-1> shp;
         std::copy(shape.begin(), shape.end() - 1, shp.begin());
         types::ndarray<T,N-1> prody(shp, __builtin__::None);
         std::transform(array.begin(), array.end(), prody.begin(), [=](types::ndarray<T,N-1> const& other) {return prod(other, axis-1);});
         return prody;
     }
 }

int main (int, const char **)
{
  typedef float                                           ScalarType;    //feel free to change this to 'double' if supported by your hardware
  typedef boost::numeric::ublas::matrix<ScalarType>       MatrixType;
  
  typedef viennacl::matrix<ScalarType, viennacl::row_major>    VCLMatrixType;
  
  std::size_t dim_large = 5;
  std::size_t dim_small = 3;
  
  //
  // Setup ublas objects and fill with data:
  //
  MatrixType ublas_A(dim_large, dim_large);
  MatrixType ublas_B(dim_small, dim_small);
  MatrixType ublas_C(dim_large, dim_small);
  MatrixType ublas_D(dim_small, dim_large);
  
  
  for (std::size_t i=0; i<ublas_A.size1(); ++i)
    for (std::size_t j=0; j<ublas_A.size2(); ++j)
      ublas_A(i,j) = static_cast<ScalarType>((i+1) + (j+1)*(i+1));

  for (std::size_t i=0; i<ublas_B.size1(); ++i)
    for (std::size_t j=0; j<ublas_B.size2(); ++j)
      ublas_B(i,j) = static_cast<ScalarType>((i+1) + (j+1)*(i+1));

  for (std::size_t i=0; i<ublas_C.size1(); ++i)
    for (std::size_t j=0; j<ublas_C.size2(); ++j)
      ublas_C(i,j) = static_cast<ScalarType>((j+2) + (j+1)*(i+1));

  for (std::size_t i=0; i<ublas_D.size1(); ++i)
    for (std::size_t j=0; j<ublas_D.size2(); ++j)
      ublas_D(i,j) = static_cast<ScalarType>((j+2) + (j+1)*(i+1));
  
  //
  // Extract submatrices using the ranges in ublas
  //
  boost::numeric::ublas::range ublas_r1(0, dim_small);                      //the first 'dim_small' entries
  boost::numeric::ublas::range ublas_r2(dim_large - dim_small, dim_large);  //the last 'dim_small' entries
  boost::numeric::ublas::matrix_range<MatrixType> ublas_A_sub1(ublas_A, ublas_r1, ublas_r1); //upper left part of A
  boost::numeric::ublas::matrix_range<MatrixType> ublas_A_sub2(ublas_A, ublas_r2, ublas_r2); //lower right part of A

  boost::numeric::ublas::matrix_range<MatrixType> ublas_C_sub(ublas_C, ublas_r1, ublas_r1); //upper left part of C
  boost::numeric::ublas::matrix_range<MatrixType> ublas_D_sub(ublas_D, ublas_r1, ublas_r1); //upper left part of D

  //
  // Setup ViennaCL objects
  //
  VCLMatrixType vcl_A(dim_large, dim_large);
  VCLMatrixType vcl_B(dim_small, dim_small);
  VCLMatrixType vcl_C(dim_large, dim_small);
  VCLMatrixType vcl_D(dim_small, dim_large);
  
  viennacl::copy(ublas_A, vcl_A);
  viennacl::copy(ublas_B, vcl_B);
  viennacl::copy(ublas_C, vcl_C);
  viennacl::copy(ublas_D, vcl_D);
  
  //
  // Extract submatrices using the ranges in ViennaCL
  //
  viennacl::range vcl_r1(0, dim_small);   //the first 'dim_small' entries
  viennacl::range vcl_r2(dim_large - dim_small, dim_large); //the last 'dim_small' entries
  viennacl::matrix_range<VCLMatrixType>   vcl_A_sub1(vcl_A, vcl_r1, vcl_r1); //upper left part of A
  viennacl::matrix_range<VCLMatrixType>   vcl_A_sub2(vcl_A, vcl_r2, vcl_r2); //lower right part of A
  
  viennacl::matrix_range<VCLMatrixType>   vcl_C_sub(vcl_C, vcl_r1, vcl_r1); //upper left part of C
  viennacl::matrix_range<VCLMatrixType>   vcl_D_sub(vcl_D, vcl_r1, vcl_r1); //upper left part of D

  //
  // Copy from ublas to submatrices and back:
  //
  
  ublas_A_sub1 = ublas_B;
  viennacl::copy(ublas_B, vcl_A_sub1);
  viennacl::copy(vcl_A_sub1, ublas_B);
  
  //
  // Addition:
  //
  
  // range to range:
  ublas_A_sub2 += ublas_A_sub2;
  vcl_A_sub2 += vcl_A_sub2;

  // range to matrix:
  ublas_B += ublas_A_sub2;
  vcl_B += vcl_A_sub2;

  
  //
  // use matrix range with matrix-matrix product:
  //
  ublas_A_sub1 += prod(ublas_C_sub, ublas_D_sub);
  vcl_A_sub1 += viennacl::linalg::prod(vcl_C_sub, vcl_D_sub);

  //
  // Print result matrices:
  //
//.........这里部分代码省略.........