本文整理汇总了C++中Tensor类的典型用法代码示例。如果您正苦于以下问题:C++ Tensor类的具体用法?C++ Tensor怎么用?C++ Tensor使用的例子?那么恭喜您, 这里精选的类代码示例或许可以为您提供帮助。
在下文中一共展示了Tensor类的15个代码示例,这些例子默认根据受欢迎程度排序。您可以为喜欢或者感觉有用的代码点赞,您的评价将有助于系统推荐出更棒的C++代码示例。
示例1: libmesh_assert
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]));
}
}
}
}
示例2: libmesh_assert
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;
}
示例3: CreateLinearSystem
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));
}
}
}
示例4: dist
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
//.........这里部分代码省略.........
示例5: embedding_bag_cpu
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);
}
);
}
}
示例6: MapPlan
inline void MapPlan(Tensor<gpu,dim> _dst, const expr::Plan<E> &plan){
cuda::MapPlan<Saver>( _dst.FlatTo2D(), plan );
}
示例7: ASSERT
void Tensor<Dtype>::ShareMem(const Tensor& other) {
ASSERT(count_ == other.count(), "");
mem_ = other.mem();
}
示例8: copyTo
void Tensor::copyTo(Tensor& dst)
{
void* p = map();
dst.reshape((const char*)p, shape_, format_);
unMap();
}
示例9: libmesh_assert
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]));
//.........这里部分代码省略.........
示例10: test_fundamentals
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;
//.........这里部分代码省略.........
示例11: TEST_CASE
#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));
}
示例12: upsample_nearest1d_cpu
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;
}
示例13: test_resize
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);
}
示例14: Vector
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
//.........这里部分代码省略.........
示例15: size
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");
};