C++ Tensor类代码示例

void HPCoarsenTest::add_projection(const System & system,
                                   const Elem * elem,
                                   unsigned int var)
{
  // If we have children, we need to add their projections instead
  if (!elem->active())
    {
      libmesh_assert(!elem->subactive());
      for (unsigned int c = 0; c != elem->n_children(); ++c)
        this->add_projection(system, elem->child(c), var);
      return;
    }

  // The DofMap for this system
  const DofMap & dof_map = system.get_dof_map();

  // The type of finite element to use for this variable
  const FEType & fe_type = dof_map.variable_type (var);

  const FEContinuity cont = fe->get_continuity();

  fe->reinit(elem);

  dof_map.dof_indices(elem, dof_indices, var);

  const unsigned int n_dofs =
    cast_int<unsigned int>(dof_indices.size());

  FEInterface::inverse_map (system.get_mesh().mesh_dimension(),
                            fe_type, coarse, *xyz_values, coarse_qpoints);

  fe_coarse->reinit(coarse, &coarse_qpoints);

  const unsigned int n_coarse_dofs =
    cast_int<unsigned int>(phi_coarse->size());

  if (Uc.size() == 0)
    {
      Ke.resize(n_coarse_dofs, n_coarse_dofs);
      Ke.zero();
      Fe.resize(n_coarse_dofs);
      Fe.zero();
      Uc.resize(n_coarse_dofs);
      Uc.zero();
    }
  libmesh_assert_equal_to (Uc.size(), phi_coarse->size());

  // Loop over the quadrature points
  for (unsigned int qp=0; qp<qrule->n_points(); qp++)
    {
      // The solution value at the quadrature point
      Number val = libMesh::zero;
      Gradient grad;
      Tensor hess;

      for (unsigned int i=0; i != n_dofs; i++)
        {
          dof_id_type dof_num = dof_indices[i];
          val += (*phi)[i][qp] *
            system.current_solution(dof_num);
          if (cont == C_ZERO || cont == C_ONE)
            grad.add_scaled((*dphi)[i][qp],system.current_solution(dof_num));
          // grad += (*dphi)[i][qp] *
          //  system.current_solution(dof_num);
          if (cont == C_ONE)
            hess.add_scaled((*d2phi)[i][qp], system.current_solution(dof_num));
          // hess += (*d2phi)[i][qp] *
          //  system.current_solution(dof_num);
        }

      // The projection matrix and vector
      for (unsigned int i=0; i != Fe.size(); ++i)
        {
          Fe(i) += (*JxW)[qp] *
            (*phi_coarse)[i][qp]*val;
          if (cont == C_ZERO || cont == C_ONE)
            Fe(i) += (*JxW)[qp] *
              (grad*(*dphi_coarse)[i][qp]);
          if (cont == C_ONE)
            Fe(i) += (*JxW)[qp] *
              hess.contract((*d2phi_coarse)[i][qp]);
          // Fe(i) += (*JxW)[qp] *
          //  (*d2phi_coarse)[i][qp].contract(hess);

          for (unsigned int j=0; j != Fe.size(); ++j)
            {
              Ke(i,j) += (*JxW)[qp] *
                (*phi_coarse)[i][qp]*(*phi_coarse)[j][qp];
              if (cont == C_ZERO || cont == C_ONE)
                Ke(i,j) += (*JxW)[qp] *
                  (*dphi_coarse)[i][qp]*(*dphi_coarse)[j][qp];
              if (cont == C_ONE)
                Ke(i,j) += (*JxW)[qp] *
                  ((*d2phi_coarse)[i][qp].contract((*d2phi_coarse)[j][qp]));
            }
        }
    }
}

void MeshFunction::hessian (const Point& p,
                            const Real,
                            std::vector<Tensor>& output,
                            const std::set<subdomain_id_type>* subdomain_ids)
{
  libmesh_assert (this->initialized());

  const Elem* element = this->find_element(p,subdomain_ids);

  if (!element)
    {
      output.resize(0);
    }
  else
    {
      // resize the output vector to the number of output values
      // that the user told us
      output.resize (this->_system_vars.size());


      {
        const unsigned int dim = element->dim();


        /*
         * Get local coordinates to feed these into compute_data().
         * Note that the fe_type can safely be used from the 0-variable,
         * since the inverse mapping is the same for all FEFamilies
         */
        const Point mapped_point (FEInterface::inverse_map (dim,
                                                            this->_dof_map.variable_type(0),
                                                            element,
                                                            p));

        std::vector<Point> point_list (1, mapped_point);

        // loop over all vars
        for (unsigned int index=0; index < this->_system_vars.size(); index++)
          {
            /*
             * the data for this variable
             */
            const unsigned int var = _system_vars[index];
            const FEType& fe_type = this->_dof_map.variable_type(var);

            UniquePtr<FEBase> point_fe (FEBase::build(dim, fe_type));
            const std::vector<std::vector<RealTensor> >& d2phi =
              point_fe->get_d2phi();
            point_fe->reinit(element, &point_list);

            // where the solution values for the var-th variable are stored
            std::vector<dof_id_type> dof_indices;
            this->_dof_map.dof_indices (element, dof_indices, var);

            // interpolate the solution
            Tensor hess;

            for (unsigned int i=0; i<dof_indices.size(); i++)
              hess.add_scaled(d2phi[i][0], this->_vector(dof_indices[i]));

            output[index] = hess;
          }
      }
    }

  // all done
  return;
}

void Tensor::CreateLinearSystem(vector<double>& B_vec, Matrix& A_matrix, 
						Tensor& X, Tensor& A, Tensor& B,
						vector<int>& mult_modesX, vector<int>& mult_modesA)
{
	// fake multiply x and A together to create B, creating the linear system in the process

	assert(mult_modesX.size() == mult_modesA.size());

	if (X.Order() == mult_modesX.size() && A.Order() == mult_modesA.size())
	{
		assert(0);
	}

	int numMultElements = 1;	
	vector<int> mult_dims(mult_modesX.size(), 0);

	for (int i = 0; i < mult_modesX.size(); ++i) 
	{
		assert(X.Dim(mult_modesX[i]) == A.Dim(mult_modesA[i]));
		mult_dims[i] = X.Dim(mult_modesX[i]);
		numMultElements = numMultElements * mult_dims[i];
	}
	vector<int> mult_offsets;
	ComputeOffsets(mult_offsets, mult_dims);
	int result_order = X.Order() + A.Order() - mult_modesX.size() - mult_modesA.size(); 

	if (result_order == 0)
		assert(0);

	vector<int> result_dims;

	vector<int> free_modesX;
	vector<int> free_modesA;


	// find free indices from X
	for (int i = 0; i < X.Order(); ++i)
	{
		if (!VectorPlus::Contains(mult_modesX, i))
		{
			free_modesX.push_back(i);
		}
	}

	// find free indices from A
	for (int i = 0; i < A.Order(); ++i)
	{
		if (!VectorPlus::Contains(mult_modesA, i))
		{
			free_modesA.push_back(i);
		}
	}
	vector<int> a_mat_dims = VectorPlus::CreatePair(B.NumElements(), X.NumElements());
	A_matrix.Initialize(a_mat_dims);
	B_vec.reserve(B.NumElements());
	// fill in elements from result tensor

	FastIndexer B_indexer(B.Dims());

	for (int n = 0; n < B.NumElements(); ++n)
	{
		B_vec.push_back(B.At(n));

		vector<int>& indices = B_indexer.GetNext();
		vector<int> free_indicesX;
		vector<int> free_indicesA;
	//	B.ComputeIndexArray(indices, n);

		for (int i = 0; i < B.Order(); ++i)
		{
			if (!VectorPlus::Contains(mult_modesX, i))
				free_indicesX.push_back(indices[i]);
			else
				free_indicesA.push_back(indices[i]);
		}

		// sum over elementwise products of mult-mode elements
		double temp_sum = 0;
		FastIndexer mult_indexer(mult_dims);
		for (int k = 0; k < numMultElements; ++k)
		{
			vector<int>& mult_indices = mult_indexer.GetNext();
		//	ComputeIndexArray(mult_indices, mult_offsets, k);

			vector<int> indicesX; 
			vector<int> indicesA;

			MergeIndices(indicesX, mult_modesX, free_modesX, mult_indices, free_indicesX);
			MergeIndices(indicesA, mult_modesA, free_modesA, mult_indices, free_indicesA);

			
			A_matrix.Set(n, X.ComputeIndex(indicesX), A.At(indicesA));
		}
	}
}

void TrigonometricPathVessel::finish( const std::vector<double>& buffer ) {
  // Store the data calculated during mpi loop
  StoreDataVessel::finish( buffer );
  // Get current value of all arguments
  for(unsigned i=0; i<cargs.size(); ++i) cargs[i]=mymap->getArgument(i);

  // Determine closest and second closest point to current position
  double lambda=mymap->getLambda();
  std::vector<double> dist( getNumberOfComponents() ), dist2( getNumberOfComponents() );;
  retrieveSequentialValue( 0, false, dist );
  retrieveSequentialValue( 1, false, dist2 );
  iclose1=getStoreIndex(0); iclose2=getStoreIndex(1);
  double mindist1=dist[0], mindist2=dist2[0];
  if( lambda>0.0 ) {
    mindist1=-std::log( dist[0] ) / lambda;
    mindist2=-std::log( dist2[0] ) / lambda;
  }
  if( mindist2<mindist1 ) {
    double tmp=mindist1; mindist1=mindist2; mindist2=tmp;
    iclose1=getStoreIndex(1); iclose2=getStoreIndex(0);
  }
  for(unsigned i=2; i<getNumberOfStoredValues(); ++i) {
    retrieveSequentialValue( i, false, dist );
    double ndist=dist[0];
    if( lambda>0.0 ) ndist=-std::log( dist[0] ) / lambda;
    if( ndist<mindist1 ) {
      mindist2=mindist1; iclose2=iclose1;
      mindist1=ndist; iclose1=getStoreIndex(i);
    } else if( ndist<mindist2 ) {
      mindist2=ndist; iclose2=getStoreIndex(i);
    }
  }
  // And find third closest point
  int isign = iclose1 - iclose2;
  if( isign>1 ) isign=1; else if( isign<-1 ) isign=-1;
  int iclose3 = iclose1 + isign; double v2v2;
  // We now have to compute vectors connecting the three closest points to the
  // new point
  double v1v1 = (mymap->getReferenceConfiguration( iclose1 ))->calculate( mymap->getPositions(), mymap->getPbc(), mymap->getArguments(), mypack1, true );
  double v3v3 = (mymap->getReferenceConfiguration( iclose2 ))->calculate( mymap->getPositions(), mymap->getPbc(), mymap->getArguments(), mypack3, true );
  if( iclose3<0 || iclose3>=mymap->getFullNumberOfTasks() ) {
    ReferenceConfiguration* conf2=mymap->getReferenceConfiguration( iclose1 );
    v2v2=(mymap->getReferenceConfiguration( iclose2 ))->calc( conf2->getReferencePositions(), mymap->getPbc(), mymap->getArguments(),
         conf2->getReferenceArguments(), mypack2, true );
    (mymap->getReferenceConfiguration( iclose2 ))->extractDisplacementVector( conf2->getReferencePositions(), mymap->getArguments(),
        conf2->getReferenceArguments(), false, projdir );
  } else {
    ReferenceConfiguration* conf2=mymap->getReferenceConfiguration( iclose3 );
    v2v2=(mymap->getReferenceConfiguration( iclose1 ))->calc( conf2->getReferencePositions(), mymap->getPbc(), mymap->getArguments(),
         conf2->getReferenceArguments(), mypack2, true );
    (mymap->getReferenceConfiguration( iclose1 ))->extractDisplacementVector( conf2->getReferencePositions(), mymap->getArguments(),
        conf2->getReferenceArguments(), false, projdir );
  }

  // Stash derivatives of v1v1
  for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) mypack1_stashd_args[i]=mypack1.getArgumentDerivative(i);
  if( mymap->getNumberOfAtoms()>0 ) {
    ReferenceAtoms* at = dynamic_cast<ReferenceAtoms*>( mymap->getReferenceConfiguration( iclose1 ) );
    const std::vector<double> & displace( at->getDisplace() );
    for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) {
      mypack1_stashd_atoms[i]=mypack1.getAtomDerivative(i); mypack1.getAtomsDisplacementVector()[i] /= displace[i];
    }
  }
  // Calculate the dot product of v1 with v2
  double v1v2 = (mymap->getReferenceConfiguration(iclose1))->projectDisplacementOnVector( projdir, mymap->getArguments(), cargs, mypack1 );

  // This computes s value
  double spacing = mymap->getPropertyValue( iclose1, (mymap->property.begin())->first ) - mymap->getPropertyValue( iclose2, (mymap->property.begin())->first );
  double root = sqrt( v1v2*v1v2 - v2v2 * ( v1v1 - v3v3) );
  dx = 0.5 * ( (root + v1v2) / v2v2 - 1.);
  double path_s = mymap->getPropertyValue(iclose1, (mymap->property.begin())->first ) + spacing * dx; sp->set( path_s );
  double fact = 0.25*spacing / v2v2;
  // Derivative of s wrt arguments
  for(unsigned i=0; i<mymap->getNumberOfArguments(); ++i) {
    sp->setDerivative( i, fact*( mypack2.getArgumentDerivative(i) + (v2v2 * (-mypack1_stashd_args[i] + mypack3.getArgumentDerivative(i))
                                 + v1v2*mypack2.getArgumentDerivative(i) )/root ) );
  }
  // Derivative of s wrt atoms
  unsigned narg=mymap->getNumberOfArguments(); Tensor vir; vir.zero(); fact = 0.5*spacing / v2v2;
  if( mymap->getNumberOfAtoms()>0 ) {
    for(unsigned i=0; i<mymap->getNumberOfAtoms(); ++i) {
      Vector ader = fact*(( v1v2*mypack1.getAtomDerivative(i) + 0.5*v2v2*(-mypack1_stashd_atoms[i] + mypack3.getAtomDerivative(i) ) )/root + mypack1.getAtomDerivative(i) );
      for(unsigned k=0; k<3; ++k) sp->setDerivative( narg+3*i+k, ader[k] );
      vir-=Tensor( mymap->getPosition(i), ader );
    }
    // Set the virial
    unsigned nbase=narg+3*mymap->getNumberOfAtoms();
    for(unsigned i=0; i<3; ++i) for(unsigned j=0; j<3; ++j) sp->setDerivative( nbase+3*i+j, vir(i,j) );
  }
  // Now compute z value
  ReferenceConfiguration* conf2=mymap->getReferenceConfiguration( iclose1 );
  double v4v4=(mymap->getReferenceConfiguration( iclose2 ))->calc( conf2->getReferencePositions(), mymap->getPbc(), mymap->getArguments(),
              conf2->getReferenceArguments(), mypack2, true );
  // Extract vector connecting frames
  (mymap->getReferenceConfiguration( iclose2 ))->extractDisplacementVector( conf2->getReferencePositions(), mymap->getArguments(),
      conf2->getReferenceArguments(), false, projdir );
  // Calculate projection of vector on line connnecting frames
  double proj = (mymap->getReferenceConfiguration(iclose1))->projectDisplacementOnVector( projdir, mymap->getArguments(), cargs, mypack1 );
  double path_z = v1v1 + dx*dx*v4v4 - 2*dx*proj;
  // Derivatives for z path
//.........这里部分代码省略.........

std::tuple<Tensor, Tensor, Tensor, Tensor>
embedding_bag_cpu(const Tensor &weight, const Tensor &indices__,
                  const Tensor &offsets__, const bool scale_grad_by_freq,
                  const int64_t mode, bool sparse) {
  auto indices_arg = TensorArg(indices__, "indices__", 1);
  checkScalarType("embedding_bag", indices_arg, kLong);
  auto offsets_arg = TensorArg(offsets__, "offsets__", 1);
  checkScalarType("embedding_bag", offsets_arg, kLong);
  Tensor indices = indices__.contiguous();
  Tensor offsets = offsets__.contiguous();
  auto weight_arg = TensorArg(weight, "weight", 1);
  checkScalarTypes("embedding_bag", weight_arg, {kFloat, kDouble});

  auto bag_size = at::zeros(offsets.sizes(), indices.type());
  make_bag_size(offsets, indices, mode, bag_size);

  // If the last entries are empty, that the last offsets are irrelevant as they
  // won't change anything in the assignment of ID -> bag, but index_add would
  // throw out of bounds error. So to keep it simple we just add one more
  // entry to the end then get rid of it after make_offset2bag.
  auto offset2bag = at::zeros(
     {indices.sizes()[0] + 1}, indices__.type()); // offset2bag = [0 0 0 0 0]

  make_offset2bag(offsets, indices, offset2bag);

  offset2bag.resize_({indices.sizes()[0]});

  auto output = at::zeros({offsets.size(0), weight.size(1)}, weight.type());

  if (mode == MODE_MEAN || mode == MODE_SUM) {
    if (weight.type().scalarType() == kFloat) {
      index_select_add<float>(indices, offset2bag, weight, output);
    } else if (weight.type().scalarType() == kDouble) {
      index_select_add<double>(indices, offset2bag, weight, output);
    }
    auto ret = apply_bag_size(offsets, indices, mode, output, bag_size);
    return std::tuple<Tensor, Tensor, Tensor, Tensor>(ret, offset2bag, bag_size, bag_size);
  } else { // MODE_MAX
    return AT_DISPATCH_FLOATING_TYPES_AND_HALF(
      weight.type(), "embedding_bag_cpu_max", [&]() {
        return embedding_bag_cpu_max<scalar_t>(weight, indices, offset2bag, output, bag_size, offsets);
      }
    );
  }
}

inline void MapPlan(Tensor<gpu,dim> _dst, const expr::Plan<E> &plan){
  cuda::MapPlan<Saver>( _dst.FlatTo2D(), plan );
}

void Tensor<Dtype>::ShareMem(const Tensor& other) {
  ASSERT(count_ == other.count(), "");
  mem_ = other.mem();
}

void Tensor::copyTo(Tensor& dst)
{
    void* p = map();
    dst.reshape((const char*)p, shape_, format_);
    unMap();
}

void MeshFunction::hessian (const Point& p,
                            const Real,
                            std::vector<Tensor>& output)
{
  libmesh_assert (this->initialized());

  /* Ensure that in the case of a master mesh function, the
     out-of-mesh mode is enabled either for both or for none.  This is
     important because the out-of-mesh mode is also communicated to
     the point locator.  Since this is time consuming, enable it only
     in debug mode.  */
#ifdef DEBUG
  if (this->_master != NULL)
    {
      const MeshFunction* master =
        cast_ptr<const MeshFunction*>(this->_master);
      if(_out_of_mesh_mode!=master->_out_of_mesh_mode)
        libmesh_error_msg("ERROR: If you use out-of-mesh-mode in connection with master mesh " \
                          << "functions, you must enable out-of-mesh mode for both the master and the slave mesh function.");
    }
#endif

  // locate the point in the other mesh
  const Elem* element = this->_point_locator->operator()(p);

  // If we have an element, but it's not a local element, then we
  // either need to have a serialized vector or we need to find a
  // local element sharing the same point.
  if (element &&
      (element->processor_id() != this->processor_id()) &&
      _vector.type() != SERIAL)
    {
      // look for a local element containing the point
      std::set<const Elem*> point_neighbors;
      element->find_point_neighbors(p, point_neighbors);
      element = NULL;
      std::set<const Elem*>::const_iterator       it  = point_neighbors.begin();
      const std::set<const Elem*>::const_iterator end = point_neighbors.end();
      for (; it != end; ++it)
        {
          const Elem* elem = *it;
          if (elem->processor_id() == this->processor_id())
            {
              element = elem;
              break;
            }
        }
    }

  if (!element)
    {
      output.resize(0);
    }
  else
    {
      // resize the output vector to the number of output values
      // that the user told us
      output.resize (this->_system_vars.size());


      {
        const unsigned int dim = this->_eqn_systems.get_mesh().mesh_dimension();


        /*
         * Get local coordinates to feed these into compute_data().
         * Note that the fe_type can safely be used from the 0-variable,
         * since the inverse mapping is the same for all FEFamilies
         */
        const Point mapped_point (FEInterface::inverse_map (dim,
                                                            this->_dof_map.variable_type(0),
                                                            element,
                                                            p));

        std::vector<Point> point_list (1, mapped_point);

        // loop over all vars
        for (unsigned int index=0; index < this->_system_vars.size(); index++)
          {
            /*
             * the data for this variable
             */
            const unsigned int var = _system_vars[index];
            const FEType& fe_type = this->_dof_map.variable_type(var);

            AutoPtr<FEBase> point_fe (FEBase::build(dim, fe_type));
            const std::vector<std::vector<RealTensor> >& d2phi =
              point_fe->get_d2phi();
            point_fe->reinit(element, &point_list);

            // where the solution values for the var-th variable are stored
            std::vector<dof_id_type> dof_indices;
            this->_dof_map.dof_indices (element, dof_indices, var);

            // interpolate the solution
            Tensor hess;

            for (unsigned int i=0; i<dof_indices.size(); i++)
              hess.add_scaled(d2phi[i][0], this->_vector(dof_indices[i]));

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

bool
test_fundamentals(Index const dimension)
{
  bool
  passed = true;

  Index const
  number_components = integer_power(dimension, Tensor::ORDER);

  std::vector<Scalar> const
  X = generate_sequence<Scalar>(number_components, 1.0, 1.0);

  // Test constructor with pointer
  Tensor const
  A(dimension, &X[0]);

  // Test copy constructor
  Tensor
  B = A;

  Tensor
  C;

  // Test copy assignment
  C = B - A;

  Scalar
  error = norm_f(C);

  bool const
  copy_assigned = error <= machine_epsilon<Scalar>();
  passed = passed && copy_assigned;

  // Test fill with pointer
  B.fill(&X[0]);

  C = B - A;

  error = norm_f(C);

  bool const
  filled_pointer = error <= machine_epsilon<Scalar>();
  passed = passed && filled_pointer;

  std::vector<Scalar> const
  Y = generate_sequence<Scalar>(number_components, -1.0, -1.0);

  C.fill(&Y[0]);

  // Test increment
  C += A;

  error = norm_f(C);

  bool const
  incremented = error <= machine_epsilon<Scalar>();
  passed = passed && incremented;

  C.fill(&X[0]);

  // Test decrement
  C -= A;

  error = norm_f(C);

  bool const
  decremented = error <= machine_epsilon<Scalar>();
  passed = passed && decremented;

#ifdef HAVE_INTREPID_KOKKOSCORE
  //test Tensor fill and create for Kokkos data types
  Kokkos::View<Scalar *, Kokkos::DefaultExecutionSpace>
  X1("X1_kokkos", dimension);

  Kokkos::View<Scalar **, Kokkos::DefaultExecutionSpace>
  X2("X2_kokkos", dimension, dimension);

  Kokkos::View<Scalar ***, Kokkos::DefaultExecutionSpace>
  X3("X3_kokkos", dimension, dimension, dimension);

  Kokkos::View<Scalar ****, Kokkos::DefaultExecutionSpace>
  X4("X4_kokkos", dimension, dimension, dimension, dimension);

  Kokkos::deep_copy(X1, 3.1);
  Kokkos::deep_copy(X2, 3.2);
  Kokkos::deep_copy(X3, 3.3);
  Kokkos::deep_copy(X4, 3.4);

  Tensor
  Z(dimension); //(X1_k,0);

  Index
  rank = 0;

  Index
  temp = number_components;

  while (temp != 1) {
    temp = temp / dimension;
    rank = rank + 1;
//.........这里部分代码省略.........

#define CATCH_CONFIG_MAIN
#include "catch.hpp"

#include "ATen/ATen.h"
#include "ATen/DLConvertor.h"

#include <iostream>
#include <string.h>
#include <sstream>
#include "test_seed.h"

using namespace at;

TEST_CASE( "parallel", "[cpu]" ) {

  manual_seed(123, at::Backend::CPU);
  set_num_threads(1);

  Tensor a = rand(CPU(at::kFloat), {1,3});
  a[0][0] = 1;
  a[0][1] = 0;
  a[0][2] = 0;
  Tensor as = rand(CPU(at::kFloat), {3});
  as[0] = 1;
  as[1] = 0;
  as[2] = 0;
  REQUIRE(a.sum(0).equal(as));
}

Tensor upsample_nearest1d_cpu(const Tensor& input, IntArrayRef output_size) {
  auto output = at::empty({0}, input.options());
  upsample_nearest1d_out_cpu_template(output, input, output_size);
  return output;
}

static void test_resize()
{
  Tensor<int, 3> epsilon;
  epsilon.resize(2,3,7);
  VERIFY_IS_EQUAL(epsilon.dimension(0), 2);
  VERIFY_IS_EQUAL(epsilon.dimension(1), 3);
  VERIFY_IS_EQUAL(epsilon.dimension(2), 7);
  VERIFY_IS_EQUAL(epsilon.dimensions().TotalSize(), 2*3*7);

  const int* old_data = epsilon.data();
  epsilon.resize(3,2,7);
  VERIFY_IS_EQUAL(epsilon.dimension(0), 3);
  VERIFY_IS_EQUAL(epsilon.dimension(1), 2);
  VERIFY_IS_EQUAL(epsilon.dimension(2), 7);
  VERIFY_IS_EQUAL(epsilon.dimensions().TotalSize(), 2*3*7);
  VERIFY_IS_EQUAL(epsilon.data(), old_data);

  epsilon.resize(3,5,7);
  VERIFY_IS_EQUAL(epsilon.dimension(0), 3);
  VERIFY_IS_EQUAL(epsilon.dimension(1), 5);
  VERIFY_IS_EQUAL(epsilon.dimension(2), 7);
  VERIFY_IS_EQUAL(epsilon.dimensions().TotalSize(), 3*5*7);
  VERIFY_IS_NOT_EQUAL(epsilon.data(), old_data);
}

void ERMSD::calcMat(const std::vector<Vector> & positions,const Pbc& pbc, std::vector<Vector4d> &mat, std::vector<TensorGeneric<4,3> > &Gderi) {

  std::vector<Vector3d> pos;
  pos.resize(3*nresidues);

  std::vector<Tensor3d> deri;
  deri.resize(nresidues*9);

  std::vector<Vector> centers;
  centers.resize(nresidues);

  unsigned idx_deri = 0;

  Tensor da_dxa = (2./3.)*Tensor::identity();
  Tensor da_dxb = -(1./3.)*Tensor::identity();
  Tensor da_dxc = -(1./3.)*Tensor::identity();

  Tensor db_dxa = -(1./3.)*Tensor::identity();
  Tensor db_dxb = (2./3.)*Tensor::identity();
  Tensor db_dxc = -(1./3.)*Tensor::identity();

  // Form factors - should this be somewhere else?

  double w = 1./3.;
  Vector form_factor = Vector(2.0,2.0,1.0/0.3);

  for(unsigned res_idx=0; res_idx<natoms/3; res_idx++) {


    const unsigned at_idx = 3*res_idx;
    //center
    for (unsigned j=0; j<3; j++) {
      centers[res_idx] += w*positions[at_idx+j];
    }

    Vector3d a = delta(centers[res_idx],positions[at_idx]);
    Vector3d b = delta(centers[res_idx],positions[at_idx+1]);
    Vector3d d = crossProduct(a,b);
    double ianorm = 1./a.modulo();
    double idnorm = 1./d.modulo();

    // X vector: COM-C2
    pos[at_idx] = a*ianorm;
    // Z versor: C2 x (COM-C4/C6)
    pos[at_idx+2] = d*idnorm;
    // Y versor: Z x Y
    pos[at_idx+1] = crossProduct(pos[at_idx+2],pos[at_idx]);

    // Derivatives ////////
    Tensor3d t1 = ianorm*(Tensor::identity()-extProduct(pos[at_idx],pos[at_idx]));
    // dv1/dxa
    deri[idx_deri] = (2./3. )*t1;
    // dv1/dxb
    deri[idx_deri+3] = -(1./3.)*t1;
    // dv1/dxc
    deri[idx_deri+6] = -(1./3.)*t1;

    Tensor dd_dxa =  VcrossTensor(a,db_dxa) -VcrossTensor(b,da_dxa);
    Tensor dd_dxb =  VcrossTensor(a,db_dxb)-VcrossTensor(b,da_dxb);
    Tensor dd_dxc =  VcrossTensor(a,db_dxc)-VcrossTensor(b,da_dxc);

    // dv3/dxa
    deri[idx_deri+2] = deriNorm(d,dd_dxa);
    // dv3/dxb
    deri[idx_deri+5] = deriNorm(d,dd_dxb);
    // dv3/dxc
    deri[idx_deri+8] = deriNorm(d,dd_dxc);

    // dv2/dxa = dv3/dxa cross v1 + v3 cross dv1/dxa
    deri[idx_deri+1] =  (VcrossTensor(deri[idx_deri+2],pos[at_idx]) + \
                         VcrossTensor(pos[at_idx+2],deri[idx_deri]));
    // dv2/dxb
    deri[idx_deri+4] =  (VcrossTensor(deri[idx_deri+5],pos[at_idx]) + \
                         VcrossTensor(pos[at_idx+2],deri[idx_deri+3]));
    // dv2/dxc
    deri[idx_deri+7] =  (VcrossTensor(deri[idx_deri+8],pos[at_idx]) + \
                         VcrossTensor(pos[at_idx+2],deri[idx_deri+6]));

    idx_deri += 9;
    // End derivatives ///////

  }


  // Initialization (unnecessary?)
  for (unsigned i1=0; i1<nresidues*nresidues; i1++) {
    for (unsigned i2=0; i2<4; i2++) {
      mat[i1][i2] = 0.0;
    }
  }

  double maxdist = cutoff/form_factor[0];
  double gamma = pi/cutoff;
  unsigned idx;
  unsigned idx1 = 0;
  // Calculate mat
  for (unsigned i=0; i<nresidues; i++) {
    for (unsigned j=0; j<nresidues; j++) {

      // skip i==j
//.........这里部分代码省略.........

void Tensor<XPU, DT>::blas_vsqrt(const Tensor<XPU, DT> &in)
{ const int N = size();  CHECK_EQ (N, in.size());
  unary_vexpr_kernel<opsqrt<DT>><<<cuda_get_blocks(N), CUDA_NUM_THREADS, 0, get_calc_stream()>>>
    (N, in.dptr, dptr);
  cuda_sync_check ("cublas_vsqrt");
};