C++ prod函数代码示例

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;
          }