I have a code snippet that looks something like the code at the end of this email. I figured a copy/paste was better for future readability rather than a gist, etc.
The basic problem is that this function takes 99% of the runtime of my application. I have ended up with an Eigen solution after converting from manual loops using STL containers. While that performance is basically fine, any improvements I can get are greatly appreciated.
When compiling with clang (OSX) or gcc (linux) , this Eigen loop is about 40% faster than the STL looping solution. Great!
What is not great is the performance in MSVC. The performance in Visual Studio 2008/2010/2012 32bit/64bit (all Release, fully optimized) is significantly worse than the "naive" solution using STL containers. "Significantly" here is between 4-6 times slower. Is there something obviously weird in my code that is killing MSVC?
It seems that enabling or disabling SSE2 doesn't have any impact, which suggests that the vectorization is not working properly. The annotated assembly is 15k lines long and 7.5 MB for this file otherwise I would have attached it.
When compiling with Visual Studio, it complains about the length of some of the templates, with errors like:
warning C4503: 'Eigen::internal::unaligned_assign_impl<IsAligned>::run' : decorated name length exceeded, name was truncated
Is that borking Visual Studio?
And/or, do you have recommendations for removing the temporaries ndteuE, B_deltaE, etc? Some are needed in more than one of the following terms so I figured that was better than a copy/paste into the following lines. Or is that the wrong idea?
Yours,
Ian
void ResidualHelmholtzGeneralizedExponential::allEigen(const long double &tau, const long double &delta, HelmholtzDerivatives &derivs) throw()
{
double log_tau = log(tau), log_delta = log(delta),
_one_over_delta_ = 1/delta, _one_over_tau_ = 1/tau; // division is much slower than multiplication, so do one division here
Eigen::Map<Eigen::ArrayXd> nE(n.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> dE(d.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> tE(t.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> cE(c.data(), elements.size());
Eigen::Map<Eigen::ArrayXi> l_intE(l_int.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> l_doubleE(l_double.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> eta1E(eta1.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> eta2E(eta2.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> epsilon1E(epsilon1.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> epsilon2E(epsilon2.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> beta1E(beta1.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> beta2E(beta2.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> gamma1E(gamma1.data(), elements.size());
Eigen::Map<Eigen::ArrayXd> gamma2E(gamma2.data(), elements.size());
// ****************************************
// The u part in exp(u) and its derivatives
// ****************************************
// Set the u part of exp(u) to zero
uE.fill(0);
du_ddeltaE.fill(0);
du_dtauE.fill(0);
d2u_ddelta2E.fill(0);
d2u_dtau2E.fill(0);
d3u_ddelta3E.fill(0);
d3u_dtau3E.fill(0);
if (delta_li_in_u){
Eigen::ArrayXd u_increment = -cE*(log_delta*l_doubleE).exp(); //pow(delta,L) -> exp(L*log(delta))
uE += u_increment;
du_ddeltaE += l_doubleE*u_increment*one_over_delta;
d2u_ddelta2E += (l_doubleE-1)*l_doubleE*u_increment*one_over_delta*one_over_delta;
d3u_ddelta3E += (l_doubleE-2)*(l_doubleE-1)*l_doubleE*u_increment*one_over_delta*one_over_delta*one_over_delta;
}
if (eta1_in_u){
uE += -eta1E*(delta-epsilon1E);
du_ddeltaE += -eta1E;
}
if (eta2_in_u){
uE += -eta2E*POW2(delta-epsilon2E);
du_ddeltaE += -2*eta2E*(delta-epsilon2E);
d2u_ddelta2E += -2*eta2E;
}
if (beta1_in_u){
uE += -beta1E*(tau-gamma1E);
du_dtauE += -beta1E;
}
if (beta2_in_u){
uE += -beta2E*POW2(tau-gamma2E);
du_dtauE += -2*beta2E*(tau-gamma2E);
d2u_dtau2E += -2*beta2E;
}
Eigen::ArrayXd ndteuE = nE*exp(tE*log_tau + dE*log_delta + uE);
Eigen::ArrayXd B_deltaE = delta*du_ddeltaE + dE;
Eigen::ArrayXd B_tauE = tau*du_dtauE + tE;
Eigen::ArrayXd B_delta2E = POW2(delta)*(d2u_ddelta2E + du_ddeltaE.square()) + 2*dE*delta*du_ddeltaE + dE*(dE-1);
Eigen::ArrayXd B_tau2E = POW2(tau)*(d2u_dtau2E + du_dtauE.square()) + 2*tE*tau*du_dtauE + tE*(tE-1);
Eigen::ArrayXd B_delta3E = POW3(delta)*d3u_ddelta3E + 3*dE*POW2(delta)*d2u_ddelta2E+3*POW3(delta)*d2u_ddelta2E*du_ddeltaE+3*dE*POW2(delta*du_ddeltaE)+3*dE*(dE-1)*delta*du_ddeltaE+dE*(dE-1)*(dE-2)+POW3(delta*du_ddeltaE);
Eigen::ArrayXd B_tau3E = POW3(tau)*d3u_dtau3E + 3*tE*POW2(tau)*d2u_dtau2E+3*POW3(tau)*d2u_dtau2E*du_dtauE+3*tE*POW2(tau*du_dtauE)+3*tE*(tE-1)*tau*du_dtauE+tE*(tE-1)*(tE-2)+POW3(tau*du_dtauE);
derivs.alphar += ndteuE.sum();
derivs.dalphar_ddelta += (ndteuE*B_deltaE).sum()*one_over_delta;
derivs.dalphar_dtau += (ndteuE*B_tauE).sum()*one_over_tau;
derivs.d2alphar_ddelta2 += (ndteuE*B_delta2E).sum()*POW2(one_over_delta);
derivs.d2alphar_dtau2 += (ndteuE*B_tau2E).sum()*POW2(one_over_tau);
derivs.d2alphar_ddelta_dtau += (ndteuE*B_deltaE*B_tauE).sum()*one_over_delta*one_over_tau;
derivs.d3alphar_ddelta3 += (ndteuE*B_delta3E).sum()*POW3(one_over_delta);
derivs.d3alphar_dtau3 += (ndteuE*B_tau3E).sum()*POW3(one_over_tau);
derivs.d3alphar_ddelta2_dtau += (ndteuE*B_delta2E*B_tauE).sum()*POW2(one_over_delta)*one_over_tau;
derivs.d3alphar_ddelta_dtau2 += (ndteuE*B_deltaE*B_tau2E).sum()*one_over_delta*POW2(one_over_tau);
return;
};