C++ C2函数代码示例

 virtual void compute(const vd::Time& t)
 {
     f2 = DSAT() * C2();
     f1 = C1() * p1;
     grad(C3) = (f2);
     grad(C1) = 0 - (f1);
     grad(DSAT) = 0;
     grad(C2) = (f1) - (f2);
 }

	bool collision(const CircleHitBox &cercleBox1, const CircleHitBox &cercleBox2)
	{
		sf::Vector2f C1(cercleBox1.p), C2(cercleBox2.p);
		float d2 = (C1.x-C2.x)*(C1.x-C2.x) + (C1.y-C2.y)*(C1.y-C2.y);
		if (d2 > (cercleBox1.rayon + cercleBox2.rayon)*(cercleBox1.rayon + cercleBox2.rayon))
			return false;
		else
			return true;
	}

double SpinAdapted::CreCreDesDes::redMatrixElement(Csf c1, vector<Csf>& ladder, const SpinBlock* b)
{
  assert( build_pattern == "(((CC)(D))(D))" );
  double element = 0.0;
  int I = get_orbs()[0]; 
  int J = get_orbs()[1];
  int K = get_orbs()[2];
  int L = get_orbs()[3];

  // Must take into account how the 4-index is built from a combination of the 2-index ops
  std::vector<SpinQuantum> quantum_ladder = get_quantum_ladder().at("(((CC)(D))(D))");
  assert( quantum_ladder.size() == 3 );

  SpinQuantum deltaQuantum12 = quantum_ladder.at(0);
  SpinQuantum deltaQuantum123 = quantum_ladder.at(1);
  SpinQuantum deltaQuantum1234 = quantum_ladder.at(2);
//FIXME components != 0
  deltaQuantum[0] = deltaQuantum1234;

  // Spin quantum data for CC
  IrrepSpace sym12 = deltaQuantum12.get_symm();
  int irrep12 = deltaQuantum12.get_symm().getirrep();
  int spin12 = deltaQuantum12.get_s().getirrep();
  // Spin quantum data for (CC)D
  IrrepSpace sym123 = deltaQuantum123.get_symm();
  int irrep123 = deltaQuantum123.get_symm().getirrep();
  int spin123= deltaQuantum123.get_s().getirrep();
  // Spin quantum data for total operator
  IrrepSpace sym1234 = deltaQuantum1234.get_symm();
  int irrep1234 = deltaQuantum1234.get_symm().getirrep();
  int spin1234 = deltaQuantum1234.get_s().getirrep();

  TensorOp C1(I, 1); 
  TensorOp C2(J, 1); 
  TensorOp D3(K,-1); 
  TensorOp D4(L,-1); 

  TensorOp CC = C1.product(C2, spin12, irrep12);
  TensorOp CCD = CC.product(D3, spin123, irrep123);
  TensorOp CCDD = CCD.product(D4, spin1234, irrep1234);

//FIXME loop over deltaQuantum components
  int j = 0;
  for (int i=0; i<ladder.size(); i++)
  {
    int index = 0; double cleb=0.0;
    if (nonZeroTensorComponent(c1, deltaQuantum[j], ladder[i], index, cleb)) {
      std::vector<double> MatElements = calcMatrixElements(c1, CCDD, ladder[i]) ;
      element = MatElements[index]/cleb;
      break;
    }
    else
      continue;
  }
  return element;
}

double SpinAdapted::StackCreCreCreCre::redMatrixElement(Csf c1, vector<Csf>& ladder, const StackSpinBlock* b)
{
  assert( build_pattern == "(((CC)(C))(C))" );
  double element = 0.0;
  int I = get_orbs()[0]; 
  int J = get_orbs()[1];
  int K = get_orbs()[2];
  int L = get_orbs()[3];
  int Slaterlength = c1.det_rep.begin()->first.size();
  vector<bool> backupSlater1(Slaterlength,0), backupSlater2(Slaterlength,0);

  // Must take into account how the 4-index is built from a combination of the 2-index ops
  std::vector<SpinQuantum> quantum_ladder = get_quantum_ladder().at("(((CC)(C))(C))");
  assert( quantum_ladder.size() == 3 );

  SpinQuantum deltaQuantum12 = quantum_ladder.at(0);
  SpinQuantum deltaQuantum123 = quantum_ladder.at(1);
  SpinQuantum deltaQuantum1234 = quantum_ladder.at(2);
  deltaQuantum[0] = deltaQuantum1234;

  // Spin quantum data for CC
  IrrepSpace sym12 = deltaQuantum12.get_symm();
  int irrep12 = deltaQuantum12.get_symm().getirrep();
  int spin12 = deltaQuantum12.get_s().getirrep();
  // Spin quantum data for (CC)C
  IrrepSpace sym123 = deltaQuantum123.get_symm();
  int irrep123 = deltaQuantum123.get_symm().getirrep();
  int spin123= deltaQuantum123.get_s().getirrep();
  // Spin quantum data for total operator
  IrrepSpace sym1234 = deltaQuantum1234.get_symm();
  int irrep1234 = deltaQuantum1234.get_symm().getirrep();
  int spin1234 = deltaQuantum1234.get_s().getirrep();

  TensorOp C1(I, 1); 
  TensorOp C2(J, 1); 
  TensorOp C3(K, 1); 
  TensorOp C4(L, 1); 

  TensorOp CC = C1.product(C2, spin12, irrep12);
  TensorOp CCC = CC.product(C3, spin123, irrep123);
  TensorOp CCCC = CCC.product(C4, spin1234, irrep1234);

  for (int i=0; i<ladder.size(); i++)
  {
    int index = 0; double cleb=0.0;
    if (nonZeroTensorComponent(c1, deltaQuantum[0], ladder[i], index, cleb)) {
      std::vector<double> MatElements = calcMatrixElements(c1, CCCC, ladder[i], backupSlater1, backupSlater2) ;
      element = MatElements[index]/cleb;
      break;
    }
    else
      continue;
  }
  return element;
}

void RatioNSDistanceTransform::processRow(const BinaryPixelType *imageRow) {
    int col;
#define N1_SETMINUS_N2_COUNT  1
#define N2_SETMINUS_N1_COUNT  5
#define N1_CAP_N2_COUNT 3
    static vect n1[N1_SETMINUS_N2_COUNT]   = {{-1, 1}};
    static vect n2[N2_SETMINUS_N1_COUNT]   = {{1, 0}, {2, 0}, {2, 1}, {1, 2}, {2, 2}};
    static vect n12[N1_CAP_N2_COUNT] = {{0, 1}, {1, 1}, {0, 2}};

    for (col = 0; col < _cols; col++) {
	if (imageRow[col] == 0)
	    dtLines[0][col + 2] = 0;
	else {
	    GrayscalePixelType val;
	    GrayscalePixelType dt;
	    int k;

	    val = GRAYSCALE_MAX;
	    for (k = 0; k < N1_SETMINUS_N2_COUNT; k++) {
		assert(n1[k].y >= 0);
		assert(n1[k].y <= 2);
		assert(col + 2 - n1[k].x >= 0);
		assert(col + 2 - n1[k].x < _cols + 3);
		val = std::min(val, dtLines[n1[k].y][col + 2 - n1[k].x]);
	    }
	    assert(C1(d.num, d.den, (int) val) == d.mbf1i(d.mbf1(val)+1)+1);
	    dt = std::min((int)_dMax, d.mbf1i(d.mbf1(val)+1)+1);

	    val = GRAYSCALE_MAX;
	    for (k = 0; k < N2_SETMINUS_N1_COUNT; k++) {
		assert(n2[k].y >= 0);
		assert(n2[k].y <= 2);
		assert(col + 2 - n2[k].x >= 0);
		assert(col + 2 - n2[k].x < _cols + 3);
		val = std::min(val, dtLines[n2[k].y][col + 2 - n2[k].x]);
	    }
	    assert(C2(d.num, d.den, (int) val) == d.mbf2i(d.mbf2(val)+1)+1);
	    dt = std::min((int) dt, d.mbf2i(d.mbf2(val)+1)+1);

	    val = GRAYSCALE_MAX;
	    for (k = 0; k < N1_CAP_N2_COUNT; k++) {
		assert(n12[k].y >= 0);
		assert(n12[k].y <= 2);
		assert(col + 2 - n12[k].x >= 0);
		assert(col + 2 - n12[k].x < _cols + 3);
		val = std::min(val, dtLines[n12[k].y][col + 2 - n12[k].x]);
	    }
	    dt = std::min((int) dt, val + 1);

	    dtLines[0][col + 2] = dt;
	}
    }
    _consumer->processRow(dtLines[0]+2);
    rotate();
}

Molecule C2H4()
{
    int nAtoms = 6;

    Eigen::Vector3d C1(0.0000000000,  0.0000000000,  1.2578920000);
    Eigen::Vector3d H1(0.0000000000,  1.7454620000,  2.3427160000);
    Eigen::Vector3d H2(0.0000000000, -1.7454620000,  2.3427160000);
    Eigen::Vector3d C2(0.0000000000,  0.0000000000, -1.2578920000);
    Eigen::Vector3d H3(0.0000000000,  1.7454620000, -2.3427160000);
    Eigen::Vector3d H4(0.0000000000, -1.7454620000, -2.3427160000);

    Eigen::MatrixXd geom(3, nAtoms);
    geom.col(0) = C1.transpose();
    geom.col(1) = H1.transpose();
    geom.col(2) = H2.transpose();
    geom.col(3) = C2.transpose();
    geom.col(4) = H3.transpose();
    geom.col(5) = H4.transpose();
    Eigen::VectorXd charges(6), masses(6);
    charges << 6.0, 1.0, 1.0, 6.0, 1.0, 1.0;
    masses  << 12.00, 1.0078250, 1.0078250, 12.0, 1.0078250, 1.0078250;

    double radiusC = (1.70 * 1.20) / convertBohrToAngstrom;
    double radiusH = (1.20 * 1.20) / convertBohrToAngstrom;
    std::vector<Atom> atoms;
    atoms.push_back( Atom("Carbon",   "C", charges(0), masses(0), radiusC, C1, 1.0) );
    atoms.push_back( Atom("Hydrogen", "H", charges(1), masses(1), radiusH, H1, 1.0) );
    atoms.push_back( Atom("Hydrogen", "H", charges(2), masses(2), radiusH, H2, 1.0) );
    atoms.push_back( Atom("Carbon",   "C", charges(3), masses(3), radiusC, C2, 1.0) );
    atoms.push_back( Atom("Hydrogen", "H", charges(4), masses(4), radiusH, H3, 1.0) );
    atoms.push_back( Atom("Hydrogen", "H", charges(5), masses(5), radiusH, H4, 1.0) );

    std::vector<Sphere> spheres;
    Sphere sph1(C1, radiusC);
    Sphere sph2(H1, radiusH);
    Sphere sph3(H2, radiusH);
    Sphere sph4(C2, radiusC);
    Sphere sph5(H3, radiusH);
    Sphere sph6(H4, radiusH);
    spheres.push_back(sph1);
    spheres.push_back(sph2);
    spheres.push_back(sph3);
    spheres.push_back(sph4);
    spheres.push_back(sph5);
    spheres.push_back(sph6);

    // D2h as generated by Oxy, Oxz, Oyz
    Symmetry pGroup = buildGroup(3, 4, 2, 1);

    return Molecule(nAtoms, charges, masses, geom, atoms, spheres, pGroup);
};

// Update the coefficients associated with the patch field
void Foam::smoluchowskiJumpTFvPatchScalarField::updateCoeffs()
{
    if (updated())
    {
        return;
    }

    const fvPatchScalarField& pmu =
        patch().lookupPatchField<volScalarField, scalar>("mu");
    const fvPatchScalarField& prho =
        patch().lookupPatchField<volScalarField, scalar>("rho");
    const fvPatchField<scalar>& ppsi =
        patch().lookupPatchField<volScalarField, scalar>("psi");
    const fvPatchVectorField& pU =
        patch().lookupPatchField<volVectorField, vector>("U");

    // Prandtl number reading consistent with rhoCentralFoam
    const dictionary& thermophysicalProperties =
        db().lookupObject<IOdictionary>("thermophysicalProperties");

    dimensionedScalar Pr
    (
        dimensionedScalar::lookupOrDefault
        (
            "Pr",
            thermophysicalProperties,
            1.0
        )
    );

    Field<scalar> C2
    (
        pmu/prho
        *sqrt(ppsi*constant::mathematical::piByTwo)
        *2.0*gamma_/Pr.value()/(gamma_ + 1.0)
        *(2.0 - accommodationCoeff_)/accommodationCoeff_
    );

    Field<scalar> aCoeff(prho.snGrad() - prho/C2);
    Field<scalar> KEbyRho(0.5*magSqr(pU));

    valueFraction() = (1.0/(1.0 + patch().deltaCoeffs()*C2));
    refValue() = Twall_;
    refGrad() = 0.0;

    mixedFvPatchScalarField::updateCoeffs();
}

int main(int argc, char **argv) {
  int m = 8000;
  int k = 3600;
  int n = 3600;
  int numsteps = 1;

  Matrix<double> A = RandomMatrix<double>(m, k);
  Matrix<double> B = RandomMatrix<double>(k, n);
  Matrix<double> C1(m, n), C2(m, n);
  
  Time([&] { MatMul(A, B, C1); }, "Classical gemm");
  Time([&] { grey433_29_234::FastMatmul(A, B, C2, numsteps); }, "Fast (4, 3, 3)");
  
  // Test for correctness.
  std::cout << "Maximum relative difference: " << MaxRelativeDiff(C1, C2) << std::endl;
  
  return 0;
}

/*!
 * @param rxn  Reaction index of the current reaction. This is used
 *             as an index into vectors which have length n_total_rxn.
 * @param k    This is a vector of integer values specifying the
 *             species indices. The length of this vector species
 *             the number of different species in the description.
 *             The value of the entries are the species indices.
 *             These are used as indexes into vectors which have
 *             length n_total_species.
 *  @param order This is a vector of the same length as vector k.
 *         The order is used for the routine power(), which produces
 *         a power law expression involving the species vector.
 *  @param stoich  This is used to handle fractional stoichiometric coefficients
 *                 on the product side of irreversible reactions.
 */
void StoichManagerN::add(size_t rxn, const std::vector<size_t>& k, const vector_fp& order,
                         const vector_fp& stoich) {
    //printf ("add called\n");
    if (order.size() != k.size()) {
        throw CanteraError("StoichManagerN::add()", "size of order and species arrays differ");
    }
    if (stoich.size() != k.size()) {
        throw CanteraError("StoichManagerN::add()", "size of stoich and species arrays differ");
    }
    bool frac = false;
    for (size_t n = 0; n < stoich.size(); n++) {
        if (fmod(stoich[n], 1.0) || fmod(order[n], 1.0)) {
            frac = true;
            break;
        }
    }
    if (frac || k.size() > 3) {
        m_cn_list.push_back(C_AnyN(rxn, k, order, stoich));
    } else {
        // Try to express the reaction with unity stoichiometric
        // coefficients (by repeating species when necessary) so that the
        // simpler 'multiply' function can be used to compute the rate
        // instead of 'power'.
        std::vector<size_t> kRep;
        for (size_t n = 0; n < k.size(); n++) {
            for (size_t i = 0; i < stoich[n]; i++)
                kRep.push_back(k[n]);
        }

        switch (kRep.size()) {
        case 1:
            m_c1_list.push_back(C1(rxn, kRep[0]));
            break;
        case 2:
            m_c2_list.push_back(C2(rxn, kRep[0], kRep[1]));
            break;
        case 3:
            m_c3_list.push_back(C3(rxn, kRep[0], kRep[1], kRep[2]));
            break;
        default:
            m_cn_list.push_back(C_AnyN(rxn, k, order, stoich));
        }
    }
}

int main(int argc, char **argv) {
  int m = 2000;
  int k = 1200;
  int n = 2000;
  int numsteps = 1;

  Matrix<double> A = RandomMatrix<double>(m, k);
  Matrix<double> B = RandomMatrix<double>(k, n);
  Matrix<double> C1(m, n), C2(m, n);

  Time([&] { MatMul(A, B, C1); }, "Classical gemm");
  Time([&] { classical222_8_24::FastMatmul(A, B, C2, numsteps); },
       "Classical recursive (2, 2, 2)");

  // Test for correctness.
  std::cout << "Maximum relative difference: " << MaxRelativeDiff(C1, C2) << std::endl;

  return 0;
}

int main()
{
    Array<int,2> A(2,3), B(2,3);

    A = 0, 3, 5, 
        1, 6, 9;
    B = 2, 5, 1,
        9, 3, 4;

    Array<int,2> C(A+2*B);
    Array<int,2> C2(2,3);
    C2 =  4, 13,  7,
         19, 12, 17;
    BZTEST(count(C2 == C) == 6); 

    beginCheckAssert();
    Array<int,2> D(i*10+j);
    endCheckAssert();
}

//======================================================================
tensor EightNode_Brick_u_p::getDampingTensorS( )
{
    int C2_dim[] = {Num_Nodes, Num_Nodes};
    tensor C2(2,C2_dim,0.0);
    tensor Cc2(2,C2_dim,0.0);

    double r  = 0.0;
    double rw = 0.0;
    double s  = 0.0;
    double sw = 0.0;
    double t  = 0.0;
    double tw = 0.0;
    double weight = 0.0;
    double det_of_Jacobian = 1.0;

    tensor Jacobian;
    tensor h;
    double Qqinv = (alpha-nf)/ks + nf/kf;

    int GP_c_r, GP_c_s, GP_c_t;

    for( GP_c_r = 0 ; GP_c_r < Num_IntegrationPts; GP_c_r++ ) {
      r = pts[GP_c_r];
      rw = wts[GP_c_r];
      for( GP_c_s = 0 ; GP_c_s < Num_IntegrationPts; GP_c_s++ ) {
        s = pts[GP_c_s];
        sw = wts[GP_c_s];
        for( GP_c_t = 0 ; GP_c_t < Num_IntegrationPts; GP_c_t++ ) {
          t = pts[GP_c_t];
          tw = wts[GP_c_t];
  	      h = shapeFunction(r,s,t);
	      Jacobian = this->Jacobian_3D(r,s,t);
	      det_of_Jacobian  = Jacobian.determinant();
	      weight = rw * sw * tw * det_of_Jacobian;
	      Cc2 = h("a")*h("k");
	      C2 += Cc2*(weight*Qqinv);
	   }
	}
    }

    return C2;
}

Z H2(xpn){
  A z;
  ND2 I1;
  {I ar=a->r,an=a->n;XW;
   I bn=b0(a->p,an),wl=wr?*wd:1;
   aw=0;
   Q(bn<0,9)
   Q(ar>1,7)
   if(!wr) wl=wr=1;
   if(wn==1) aw=2;
   else Q(wl!=bn,8)
   if(wr==1&&wt!=Et){
     W(gv(wt,an))
     C2((I(*)())(!wt?(I(*)())x0:wt==Ft?(I(*)())x1:(I(*)())x2))
   }
   v=tr(wr-1,wd+1);
   u=wn;
   W(ga(t=wt,wr,an*v,wd))*z->d=an;
   C2(x3)
}}

int main(int argc, char **argv) {
  int m = 90;
  int k = 90;
  int n = 90;
  int numsteps = 1;

  Matrix<double> A = RandomMatrix<double>(m, k);
  Matrix<double> B = RandomMatrix<double>(k, n);
  Matrix<double> C1(m, n), C2(m, n);
  MatMul(A, B, C1);

  for (int numsteps = 1; numsteps <= 2; ++numsteps) {
	double lambda = DBL_EPSILON;
	std::cout << numsteps << std::endl;
	while (lambda < 1) {
	  smirnov333_20_182_approx::FastMatmul(A, B, C2, numsteps, lambda);
	  std::cout << lambda << ", " << MaxRelativeDiff(C1, C2) << "; ";
	  lambda *= 2;
	}
	std::cout << std::endl;
  }

  return 0;
}

cString cPlainKeyVia::PrintKeyNr(void)
{
  char tmp[12];
  const char *kn=tmp;
  switch(keynr) {
    case MBC3('T','P','S'):
      kn="TPS"; break;
    case MBC3('M','K',0): case MBC3('M','K',1): case MBC3('M','K',2):
    case MBC3('M','K',3): case MBC3('M','K',4): case MBC3('M','K',5):
    case MBC3('M','K',6): case MBC3('M','K',7): case MBC3('M','K',8):
    case MBC3('M','K',9):
      snprintf(tmp,sizeof(tmp),"TPSMK%d",C3(keynr)); break;
    default:
      {
      char c2=C2(keynr);
      if(c2=='D' || c2=='P' || c2=='X' || c2=='C' || c2=='E' || c2=='T')
        snprintf(tmp,sizeof(tmp),"%c%01X",c2,keynr & 0x0f);
      else
        snprintf(tmp,sizeof(tmp),"%02X",keynr);
      break;
      }
    }
  return kn;
}