Re: [eigen] Optimization advice for a specific expression

```I am reattaching the code, as it did not get sent correctly:

```
```#include <iostream>
#include <ctime>
#include <chrono>

#include <Eigen/Core>
#include <Eigen/Geometry>

// @InProceedings{Euler2006,
//   Title                    = {Fast Computation Of Inertia Through
//                               Affinely Extended Euler And Tensor},
//   Author                   = {DiCarlo, Antonio and Paoluzzi, Alberto},
//   Year                     = {2006}}

template<class REAL, class INT>
class InertiaTensor
{

public:

// ----------------------------------------------------------------------------------------- Constructor
/// Constructor without parameters.
InertiaTensor()
: dens_(0.0)
{
Eplus << 2 , 1 , 1 , 5,
1 , 2 , 1 , 5,
1 , 1 , 2 , 5,
5 , 5 , 5 , 20; // 120 * E+
}

// ------------------------------------------------------------------------------------------- Accessors
/// Set the density of the volume.
/// \param dens         Density of the volume
inline void setDensity(const REAL dens){ dens_ = dens; };

/// Returns the mass of the volume
/// \return The volume mass, a scalar
inline REAL  getMass()  { return mass_; };

/// Returns the position of the centre of mass of the volume (defined in the same coordinates as the mesh)
/// \return A 3x1 array
inline REAL* getrG()  { return rG_; };

/// Returns the inertia tensor of the volume about its centre of mass (defined in the same coordinates as the mesh)
/// \return A 3x3 array
inline REAL* getInertiaTensor()  { return J_; };

// ------------------------------------------------------------------------------------------- Begin Vol
/// Starts the evaluation of the volumetric properties of a new volume.
/// Clears the integrals, mass, centre of mass position, and inertia tensor.
void beginVolume()
{
InertiaTensor4x4.setZero();
}

// --------------------------------------------------------------------------------------------- End Vol
/// Finalizes the volume integration and evaluates the physical properties.
void endVolume()
{
computeMassProperties();
}

/// Adds a new face to the integration process.
/// \param vertex_data  A 3xnum_vertices array with the x, y, and z mesh coordinates of the face vertices.
/// \param normal_data  A 3x1 array with the x, y, and z mesh coordinates of the face normal. May be set to NULL if the face is a triangle.
/// \param num_vertices     Number of vertices in the face.
void addFace(const REAL* vertex_data, const REAL* normal_data, const INT num_vertices)
{
// Check that we have at least 2 vertices per face
assert(num_vertices == 3);

Eigen::Map<const Eigen::Matrix<REAL, 3, 3>> v(vertex_data);
Eigen::Matrix<REAL, 4, 4> G;

// G = [v0 v1 v2 0]
//     [ 0  0  0 1]
G.col(3).setZero();
G.row(3).setZero();
G(3, 3) = 1.0;

Eigen::Matrix<REAL, 4, 4>  temp;
Eigen::Vector4d w = G.template leftCols<3>().rowwise().sum();
temp << G.template block<3, 3>(0, 0) * G.template block<3, 3>(0, 0).transpose() + w.template head<3>() * w.transpose().template head<3>(), 5 * w.template head<3>(),

// I4x4 = det(G) G E+ Gᵀ
InertiaTensor4x4.noalias() +=  temp * (G.template block<3,3>(0,0).determinant() / 120.0);

// Initial implementation
// InertiaTensor4x4.noalias() += (G.template block<3,3>(0,0).determinant() / 120.0) * G * Eplus * G.transpose();
}

private:

REAL dens_;             // Density
REAL mass_;             // Mass of the volume
REAL rG_[3];        // Position of the centre of mass
REAL J_[9];             // Inertia tensor

Eigen::Matrix<REAL, 4, 4> Eplus, InertiaTensor4x4;

// ----------------------------------------------------------------------------- Compute mass properties
void computeMassProperties()
{
// Evaluate mass
mass_ = InertiaTensor4x4(3, 3) * dens_;
}
};

template<class F>
void massProps(const char *object, const char* method, unsigned int niters, double dens, double *mass1, double *J1, double *rg1, int nf, const F &feval) {
std::cout << "--------------------------------------------------" << std::endl;
std::cout << object << ": " << nf << " faces." << std::endl;
std::cout << "--------------------------------------------------" << std::endl;

std::chrono::high_resolution_clock::time_point t1 = std::chrono::high_resolution_clock::now();
for (unsigned int iters = 0; iters < niters; iters++)
{
feval();
}
std::chrono::high_resolution_clock::time_point t2 = std::chrono::high_resolution_clock::now();
std::chrono::duration<double> time_span1 = std::chrono::duration_cast<std::chrono::duration<double>>(t2 - t1);

std::cout << object << " - " << method << std::endl;
std::cout << "Elapsed time: " << time_span1.count() << " seconds." << std::endl;
std::cout << "Total mass: " << *mass1 << std::endl;
std::cout << std::endl;
}

template <class DT>
int computeMassPropertiesIT(DT *v, int *i, int nf, const double dens,
double *mass, double *r, double *J)
{
InertiaTensor<double, int> vp;
vp.setDensity(dens);

double faceverts[9];

Eigen::Map<Eigen::Matrix<DT, 3 , Eigen::Dynamic>> verts(v, 3, *std::max_element(i, i + nf) + 1);

vp.beginVolume();

for (int face = 0; face < nf; face++)
{
Eigen::Map<Eigen::Matrix3d> points(faceverts);

points.col(0) = verts.col(i[face * 3 + 0]).template cast<double>();
points.col(1) = verts.col(i[face * 3 + 1]).template cast<double>();
points.col(2) = verts.col(i[face * 3 + 2]).template cast<double>();

}

vp.endVolume();

*mass = vp.getMass();
Eigen::Map<Eigen::Vector3d>(r, 3, 1) = Eigen::Map<Eigen::Vector3d>(vp..getrG(), 3, 1);
Eigen::Map<Eigen::Matrix3d>(J, 3, 3) = Eigen::Map<Eigen::Matrix3d>(vp..getInertiaTensor(), 3, 3);
}

int main(int argc, char *argv[])
{
// Definition of a tetrahedron ____________________________________________________________

// Number of faces
int nf = 4;

// Array of normals
double n[12] = {
1.0 / sqrt(3.0), 1.0 / sqrt(3.0), 1.0 / sqrt(3.0),
0.0, -1.0, 0.0,
0.0, 0.0, -1.0,
-1.0, 0.0, 0.0,
};

// Array of vertices (x, y, z coordinates of each vertex)
double v[12] = {
0.0, 0.0, 0.0,
1.0, 0.0, 0.0,
0.0, 1.0, 0.0,
0.0, 0.0, 1.0
};

// Array of vertex indices (vertices in each face)
int i[12] = {
1, 2, 3,
0, 1, 3,
0, 2, 1,
0, 3, 2
};

// Evaluation of mass properties __________________________________________________________

double rg1[3], J1[9];
double mass1;
const double dens = 7780.0;

const unsigned int niters_pyramid = 1000000;

std::cout << "Using Eigen version: " << EIGEN_WORLD_VERSION << "." << EIGEN_MAJOR_VERSION << "." << EIGEN_MINOR_VERSION << "\n";

massProps("PYRAMID", "DICARLO", niters_pyramid, dens, &mass1, J1, rg1, nf,
[&](){ computeMassPropertiesIT(v, i, nf, dens, &mass1, rg1, J1);});

return 0;
}
```

 Mail converted by MHonArc 2.6.19+ http://listengine.tuxfamily.org/