Skip to content
Snippets Groups Projects
Commit e39157e2 authored by Jana Pazurikova's avatar Jana Pazurikova
Browse files

Close structure - add vectorization icc pragmas

parent 253e263c
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
src/colvar/PathMSDBase.cpp
+ 8
− 0
View file @ e39157e2
......@@ -142,6 +142,7 @@ void PathMSDBase::calculate(){
// resize the list to full
if(imgVec.empty()){ // this is the signal that means: recalculate all
imgVec.resize(nframes);
#pragma simd
for(unsigned i=0;i<nframes;i++){
imgVec[i].property=indexvec[i];
imgVec[i].index=i;
......@@ -199,6 +200,7 @@ void PathMSDBase::calculate(){
for(unsigned i=rank;i<imgVec.size();i+=stride){
tmp_distances[i] = msdv[imgVec[i].index].calc_Rot(getPositions(), tmp_derivs, tmp_rotationRefClose[imgVec[i].index], true);
plumed_assert(tmp_derivs.size()==nat);
#pragma simd
for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
}
}
......@@ -207,6 +209,7 @@ void PathMSDBase::calculate(){
for(unsigned i=rank;i<imgVec.size();i+=stride){
tmp_distances[i] = msdv[imgVec[i].index].calculateWithCloseStructure(getPositions(), tmp_derivs, rotationPosClose, rotationRefClose[imgVec[i].index], drotationPosCloseDrr01, true);
plumed_assert(tmp_derivs.size()==nat);
#pragma simd
for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
if (debugClose){
double withclose = tmp_distances[i];
......@@ -226,6 +229,7 @@ void PathMSDBase::calculate(){
for(unsigned i=rank;i<imgVec.size();i+=stride){
tmp_distances[i]=msdv[imgVec[i].index].calculate(getPositions(),tmp_derivs,true);
plumed_assert(tmp_derivs.size()==nat);
#pragma simd
for(unsigned j=0;j<nat;j++) tmp_derivs2[i*nat+j]=tmp_derivs[j];
}
}
......@@ -275,12 +279,16 @@ void PathMSDBase::calculate(){
val_z_path->set(-(1./lambda)*std::log(partition));
for(unsigned j=0;j<s_path.size();j++){
// clean up
#pragma simd
for(unsigned i=0;i< derivs_s.size();i++){derivs_s[i].zero();}
// do the derivative
#pragma simd
for(const auto & it : imgVec){
double expval=it.similarity;
tmp=lambda*expval*(s_path[j]-it.property[j])/partition;
#pragma ivdep
for(unsigned i=0;i< derivs_s.size();i++){ derivs_s[i]+=tmp*it.distder[i] ;}
#pragma ivdep
if(j==0){for(unsigned i=0;i< derivs_z.size();i++){ derivs_z[i]+=it.distder[i]*expval/partition;}}
}
for(unsigned i=0;i< derivs_s.size();i++){
......
src/tools/RMSD.cpp
+ 22
− 4
View file @ e39157e2
......@@ -97,7 +97,9 @@ void RMSD::setReference(const vector<Vector> & reference){
plumed_massert(displace.empty(),"you should first clear() an RMSD object, then set a new reference");
align.resize(n,1.0/n);
displace.resize(n,1.0/n);
#pragma ivdep
for(unsigned i=0;i<n;i++) reference_center+=this->reference[i]*align[i];
#pragma ivdep
for(unsigned i=0;i<n;i++) this->reference[i]-=reference_center;
reference_center_is_calculated=true;
reference_center_is_removed=true;
......@@ -114,8 +116,10 @@ void RMSD::setAlign(const vector<double> & align, bool normalize_weights, bool r
this->align=align;
if(normalize_weights){
double w=0.0;
#pragma ivdep
for(unsigned i=0;i<n;i++) w+=this->align[i];
double inv=1.0/w;
#pragma ivdep
for(unsigned i=0;i<n;i++) this->align[i]*=inv;
}
// recalculate the center anyway
......@@ -145,8 +149,10 @@ void RMSD::setDisplace(const vector<double> & displace, bool normalize_weights){
plumed_massert(this->displace.size()==displace.size(),"mismatch in dimension of align/displace arrays");
this->displace=displace;
double w=0.0;
#pragma ivdep
for(unsigned i=0;i<n;i++) w+=this->displace[i];
double inv=1.0/w;
#pragma ivdep
if(normalize_weights){for(unsigned i=0;i<n;i++) this->displace[i]*=inv;}
}
std::vector<double> RMSD::getDisplace(){
......@@ -965,6 +971,7 @@ void RMSDCoreData::doCoreCalc(bool safe,bool alEqDis, bool only_rotation){
Vector cp; cp.zero(); if(!cpositions_is_removed)cp=cpositions;
Vector cr; cr.zero(); if(!creference_is_removed)cr=creference;
// second expensive loop: compute second moments wrt centers
#pragma ivdep
for(unsigned iat=0;iat<n;iat++){
double w=align[iat];
rr00+=dotProduct(positions[iat]-cp,positions[iat]-cp)*w;
......@@ -1075,6 +1082,7 @@ void RMSDCoreData::doCoreCalc(bool safe,bool alEqDis, bool only_rotation){
// calculate rotation matrix derivatives and components distances needed for components only when align!=displacement
if(!alEqDis)ddist_drotation.zero();
#pragma ivdep
for(unsigned iat=0;iat<n;iat++){
// components differences: this is useful externally
d[iat]=positions[iat]-cp - matmul(rotation,reference[iat]-cr);
......@@ -1105,17 +1113,20 @@ double RMSDCoreData::getDistance( bool squared){
else
dist=eigenvals[0]+rr00+rr11;
const unsigned n=static_cast<unsigned int>(reference.size());
double localDist = 0;
#pragma simd reduction(+:localDist)
for(unsigned iat=0;iat<n;iat++){
if(alEqDis){
if(safe) dist+=align[iat]*modulo2(d[iat]);
if(safe) localDist+=align[iat]*modulo2(d[iat]);
} else {
dist+=displace[iat]*modulo2(d[iat]);
localDist+=displace[iat]*modulo2(d[iat]);
}
}
if(!squared){
dist=sqrt(dist);
dist=sqrt(localDist);
distanceIsMSD=false;
}else{
dist=localDist;
distanceIsMSD=true;
}
hasDistance=true;
......@@ -1136,10 +1147,13 @@ void RMSDCoreData::doCoreCalcWithCloseStructure(bool safe,bool alEqDis, Tensor &
Tensor rotation = matmul(rotationPosClose, rotationRefClose);
#pragma simd
for (unsigned iat=0; iat<natoms; iat++){
d[iat] = positions[iat] - cp - matmul(rotation, reference[iat]-cr);
}
if (!alEqDis){
for (unsigned iat=0; iat<natoms; iat++) {
//dist = \sum w_i(x_i - cpos - R_xy * R_ay * a_i)
if (!alEqDis){
ddist_drxy += -2*displace[iat]*extProduct(matmul(d[iat], rotationRefClose), reference[iat]-cr);
}
}
......@@ -1166,6 +1180,7 @@ std::vector<Vector> RMSDCoreData::getDDistanceDPositions(){
if(!isInitialized)plumed_merror("getDPositionsDerivatives needs to initialize the coreData first!");
Vector csum;
Vector tmp1,tmp2;
#pragma ivdep
for(unsigned iat=0;iat<n;iat++){
if(alEqDis){
// there is no need for derivatives of rotation and shift here as it is by construction zero
......@@ -1184,6 +1199,7 @@ std::vector<Vector> RMSDCoreData::getDDistanceDPositions(){
}
}
#pragma ivdep
if(!alEqDis) for(unsigned iat=0;iat<n;iat++){derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcpositions-csum)*align[iat]); }
return derivatives;
......@@ -1206,6 +1222,7 @@ std::vector<Vector> RMSDCoreData::getDDistanceDReference(){
Tensor t_ddist_drr01=ddist_drr01.transpose();
// third expensive loop: derivatives
#pragma simd
for(unsigned iat=0;iat<n;iat++){
if(alEqDis){
// there is no need for derivatives of rotation and shift here as it is by construction zero
......@@ -1225,6 +1242,7 @@ std::vector<Vector> RMSDCoreData::getDDistanceDReference(){
}
}
#pragma ivdep
if(!alEqDis) for(unsigned iat=0;iat<n;iat++){derivatives[iat]= prefactor*(derivatives[iat]+(ddist_dcreference-csum)*align[iat]);}
return derivatives;
......
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment