Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Parallelization for AdjMatrices #1159

Draft
wants to merge 1 commit into
base: derivatives-from-backpropegation
Choose a base branch
from
Draft
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
144 changes: 104 additions & 40 deletions src/adjmat/AdjacencyMatrixBase.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -215,52 +215,81 @@ unsigned AdjacencyMatrixBase::retrieveNeighbours( const unsigned& current, std::

void AdjacencyMatrixBase::setupForTask( const unsigned& current, std::vector<unsigned> & indices, MultiValue& myvals ) const {
// Now retrieve bookeeping arrays
if( indices.size()!=(1+ablocks.size()+threeblocks.size()) ) indices.resize( 1+ablocks.size()+threeblocks.size() );
if( indices.size()!=(1+ablocks.size()+threeblocks.size()) ) {
indices.resize( 1+ablocks.size()+threeblocks.size() );
}

// Now get the positions
const Value* myval = getConstPntrToComponent(0);
unsigned natoms=retrieveNeighbours( current, indices ); myvals.setSplitIndex( natoms );
// const Value* myval = getConstPntrToComponent(0);
unsigned natoms=retrieveNeighbours( current, indices );
myvals.setSplitIndex( natoms );

//placing this on the stack
Vector currentAtom = ActionAtomistic::getPosition(current);
// Now retrieve everything for the third atoms
if( threeblocks.size()>0 ) {
unsigned ncells_required=0; std::vector<unsigned> cells_required( threecells.getNumberOfCells() );
threecells.addRequiredCells( threecells.findMyCell( ActionAtomistic::getPosition(current) ), ncells_required, cells_required );
unsigned ncells_required=0;
std::vector<unsigned> cells_required( threecells.getNumberOfCells() );
threecells.addRequiredCells( threecells.findMyCell( currentAtom ), ncells_required, cells_required );
threecells.retrieveAtomsInCells( ncells_required, cells_required, natoms, indices );
}
myvals.setNumberOfIndices( natoms );

// Apply periodic boundary conditions to atom positions
std::vector<std::vector<Vector> > & t_atoms( myvals.getFirstAtomDerivativeVector() ); if( t_atoms.size()!=1 ) t_atoms.resize(1);
if( t_atoms[0].size()<getNumberOfAtoms() ) t_atoms[0].resize( getNumberOfAtoms() );
for(unsigned i=0; i<natoms; ++i) t_atoms[0][i] = ActionAtomistic::getPosition(indices[i]) - ActionAtomistic::getPosition(current);
if( !nopbc ) pbcApply( t_atoms[0], natoms );
std::vector<std::vector<Vector> > & t_atoms( myvals.getFirstAtomDerivativeVector() );
//since I am usig this as a temp variable, just have to care if size is 0,
if( t_atoms.size()==0 )
t_atoms.resize(1);
if( t_atoms[0].size()<getNumberOfAtoms() )
t_atoms[0].resize( getNumberOfAtoms() );
for(unsigned i=0; i<natoms; ++i)
t_atoms[0][i] = ActionAtomistic::getPosition(indices[i]) - currentAtom;
if( !nopbc )
pbcApply( t_atoms[0], natoms );
// And collect atom position data
std::vector<Vector> & atoms( myvals.getAtomVector() );
if( atoms.size()<getNumberOfAtoms() ) atoms.resize( getNumberOfAtoms() );
for(unsigned i=0; i<natoms; ++i) atoms[ indices[i] ] = t_atoms[0][i];
if( atoms.size()<getNumberOfAtoms() )
atoms.resize( getNumberOfAtoms() );
for(unsigned i=0; i<natoms; ++i)
atoms[ indices[i] ] = t_atoms[0][i];
}

void AdjacencyMatrixBase::performTask( const std::string& controller, const unsigned& index1, const unsigned& index2, MultiValue& myvals ) const {
Vector zero; zero.zero(); plumed_dbg_assert( index2<myvals.getAtomVector().size() );
Vector zero; zero.zero();
plumed_dbg_assert( index2<myvals.getAtomVector().size() );
double weight = calculateWeight( zero, myvals.getAtomVector()[index2], myvals.getNumberOfIndices()-myvals.getSplitIndex(), myvals );
Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@gtribello
Since I do not know all the details, I have a question on calculateWeight: by git-grepping this is only called here.

And so the first argument is a fixed {0.0,0.0,0.0}, and it is only used in HBPammMatrix TopologyMatrix to do something like "pos2-0".
My suggestion is to change the signature to something like double calculateWeight(Vector , unsigned , MultiValue& ) so there is no need to copy the already calculated distance in a temporary value, moreover the signature can be that, but in the definition, you can freely add const to variables passed by value.

I am basing this suggestion also on this)

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Yes, I think if you don't pass the first argument, it would be fine. The pbc are all applied in one place because Michele Ceriotti wrote a fast algorithm to do pbc distance on many arguments at the same time, which is used here.

myvals.setValue( 0, weight );
if( fabs(weight)<epsilon ) { myvals.setValue( 0, 0 ); return; }
if( fabs(weight)<epsilon ) {
myvals.setValue( 0, 0 );
return;
}

if( !doNotCalculateDerivatives() ) {
// Update dynamic list indices for central atom
myvals.updateIndex( 0, 3*index1+0 ); myvals.updateIndex( 0, 3*index1+1 ); myvals.updateIndex( 0, 3*index1+2 );
myvals.updateIndex( 0, 3*index1+0 );
myvals.updateIndex( 0, 3*index1+1 );
myvals.updateIndex( 0, 3*index1+2 );
// Update dynamic list indices for atom forming this bond
myvals.updateIndex( 0, 3*index2+0 ); myvals.updateIndex( 0, 3*index2+1 ); myvals.updateIndex( 0, 3*index2+2 );
myvals.updateIndex( 0, 3*index2+0 );
myvals.updateIndex( 0, 3*index2+1 );
myvals.updateIndex( 0, 3*index2+2 );
// Now look after all the atoms in the third block
std::vector<unsigned> & indices( myvals.getIndices() );
for(unsigned i=myvals.getSplitIndex(); i<myvals.getNumberOfIndices(); ++i) {
myvals.updateIndex( 0, 3*indices[i]+0 ); myvals.updateIndex( 0, 3*indices[i]+1 ); myvals.updateIndex( 0, 3*indices[i]+2 );
myvals.updateIndex( 0, 3*indices[i]+0 );
myvals.updateIndex( 0, 3*indices[i]+1 );
myvals.updateIndex( 0, 3*indices[i]+2 );
}
// Update dynamic list indices for virial
unsigned base = 3*getNumberOfAtoms(); for(unsigned j=0; j<9; ++j) myvals.updateIndex( 0, base+j );
unsigned base = 3*getNumberOfAtoms();
for(unsigned j=0; j<9; ++j)
myvals.updateIndex( 0, base+j );
// And the indices for the derivatives of the row of the matrix
unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives(); std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices() );
matrix_indices[nmat_ind+0]=3*index2+0; matrix_indices[nmat_ind+1]=3*index2+1; matrix_indices[nmat_ind+2]=3*index2+2;
unsigned nmat_ind = myvals.getNumberOfMatrixRowDerivatives();
std::vector<unsigned>& matrix_indices( myvals.getMatrixRowDerivativeIndices() );
matrix_indices[nmat_ind+0]=3*index2+0;
matrix_indices[nmat_ind+1]=3*index2+1;
matrix_indices[nmat_ind+2]=3*index2+2;
myvals.setNumberOfMatrixRowDerivatives( nmat_ind+3 );
}

Expand All @@ -269,36 +298,71 @@ void AdjacencyMatrixBase::performTask( const std::string& controller, const unsi
Vector atom = myvals.getAtomVector()[index2];
myvals.setValue( 1, atom[0] ); myvals.setValue( 2, atom[1] ); myvals.setValue( 3, atom[2] );
if( !doNotCalculateDerivatives() ) {
myvals.addDerivative( 1, 3*index1+0, -1 ); myvals.addDerivative( 1, 3*index2+0, +1 );
myvals.addDerivative( 1, 3*index1+1, 0 ); myvals.addDerivative( 1, 3*index2+1, 0 );
myvals.addDerivative( 1, 3*index1+2, 0 ); myvals.addDerivative( 1, 3*index2+2, 0 );
myvals.addDerivative( 2, 3*index1+0, 0 ); myvals.addDerivative( 2, 3*index2+0, 0 );
myvals.addDerivative( 2, 3*index1+1, -1 ); myvals.addDerivative( 2, 3*index2+1, +1 );
myvals.addDerivative( 2, 3*index1+2, 0 ); myvals.addDerivative( 2, 3*index2+2, 0 );
myvals.addDerivative( 3, 3*index1+0, 0 ); myvals.addDerivative( 3, 3*index2+0, 0 );
myvals.addDerivative( 3, 3*index1+1, 0 ); myvals.addDerivative( 3, 3*index2+1, 0 );
myvals.addDerivative( 3, 3*index1+2, -1 ); myvals.addDerivative( 3, 3*index2+2, +1 );
myvals.addDerivative( 1, 3*index1+0, -1 );
myvals.addDerivative( 1, 3*index2+0, +1 );
myvals.addDerivative( 1, 3*index1+1, 0 );
myvals.addDerivative( 1, 3*index2+1, 0 );
myvals.addDerivative( 1, 3*index1+2, 0 );
myvals.addDerivative( 1, 3*index2+2, 0 );
myvals.addDerivative( 2, 3*index1+0, 0 );
myvals.addDerivative( 2, 3*index2+0, 0 );
myvals.addDerivative( 2, 3*index1+1, -1 );
myvals.addDerivative( 2, 3*index2+1, +1 );
myvals.addDerivative( 2, 3*index1+2, 0 );
myvals.addDerivative( 2, 3*index2+2, 0 );
myvals.addDerivative( 3, 3*index1+0, 0 );
myvals.addDerivative( 3, 3*index2+0, 0 );
myvals.addDerivative( 3, 3*index1+1, 0 );
myvals.addDerivative( 3, 3*index2+1, 0 );
myvals.addDerivative( 3, 3*index1+2, -1 );
myvals.addDerivative( 3, 3*index2+2, +1 );
for(unsigned k=0; k<3; ++k) {
// Update dynamic lists for central atom
myvals.updateIndex( 1, 3*index1+k ); myvals.updateIndex( 2, 3*index1+k ); myvals.updateIndex( 3, 3*index1+k );
myvals.updateIndex( 1, 3*index1+k );
myvals.updateIndex( 2, 3*index1+k );
myvals.updateIndex( 3, 3*index1+k );
// Update dynamic lists for bonded atom
myvals.updateIndex( 1, 3*index2+k ); myvals.updateIndex( 2, 3*index2+k ); myvals.updateIndex( 3, 3*index2+k );
myvals.updateIndex( 1, 3*index2+k );
myvals.updateIndex( 2, 3*index2+k );
myvals.updateIndex( 3, 3*index2+k );
}
// Add derivatives of virial
unsigned base = 3*getNumberOfAtoms();
// Virial for x
myvals.addDerivative( 1, base+0, -atom[0] ); myvals.addDerivative( 1, base+3, -atom[1] ); myvals.addDerivative( 1, base+6, -atom[2] );
myvals.addDerivative( 1, base+1, 0 ); myvals.addDerivative( 1, base+4, 0 ); myvals.addDerivative( 1, base+7, 0 );
myvals.addDerivative( 1, base+2, 0 ); myvals.addDerivative( 1, base+5, 0 ); myvals.addDerivative( 1, base+8, 0 );
myvals.addDerivative( 1, base+0, -atom[0] );
myvals.addDerivative( 1, base+3, -atom[1] );
myvals.addDerivative( 1, base+6, -atom[2] );
myvals.addDerivative( 1, base+1, 0 );
myvals.addDerivative( 1, base+4, 0 );
myvals.addDerivative( 1, base+7, 0 );
myvals.addDerivative( 1, base+2, 0 );
myvals.addDerivative( 1, base+5, 0 );
myvals.addDerivative( 1, base+8, 0 );
// Virial for y
myvals.addDerivative( 2, base+0, 0 ); myvals.addDerivative( 2, base+3, 0 ); myvals.addDerivative( 2, base+6, 0 );
myvals.addDerivative( 2, base+1, -atom[0] ); myvals.addDerivative( 2, base+4, -atom[1] ); myvals.addDerivative( 2, base+7, -atom[2] );
myvals.addDerivative( 2, base+2, 0 ); myvals.addDerivative( 2, base+5, 0 ); myvals.addDerivative( 2, base+8, 0 );
myvals.addDerivative( 2, base+0, 0 );
myvals.addDerivative( 2, base+3, 0 );
myvals.addDerivative( 2, base+6, 0 );
myvals.addDerivative( 2, base+1, -atom[0] );
myvals.addDerivative( 2, base+4, -atom[1] );
myvals.addDerivative( 2, base+7, -atom[2] );
myvals.addDerivative( 2, base+2, 0 );
myvals.addDerivative( 2, base+5, 0 );
myvals.addDerivative( 2, base+8, 0 );
// Virial for z
myvals.addDerivative( 3, base+0, 0 ); myvals.addDerivative( 3, base+3, 0 ); myvals.addDerivative( 3, base+6, 0 );
myvals.addDerivative( 3, base+1, 0 ); myvals.addDerivative( 3, base+4, 0 ); myvals.addDerivative( 3, base+7, 0 );
myvals.addDerivative( 3, base+2, -atom[0] ); myvals.addDerivative( 3, base+5, -atom[1] ); myvals.addDerivative( 3, base+8, -atom[2] );
for(unsigned k=0; k<9; ++k) { myvals.updateIndex( 1, base+k ); myvals.updateIndex( 2, base+k ); myvals.updateIndex( 3, base+k ); }
myvals.addDerivative( 3, base+0, 0 );
myvals.addDerivative( 3, base+3, 0 );
myvals.addDerivative( 3, base+6, 0 );
myvals.addDerivative( 3, base+1, 0 );
myvals.addDerivative( 3, base+4, 0 );
myvals.addDerivative( 3, base+7, 0 );
myvals.addDerivative( 3, base+2, -atom[0] );
myvals.addDerivative( 3, base+5, -atom[1] );
myvals.addDerivative( 3, base+8, -atom[2] );
for(unsigned k=0; k<9; ++k) {
myvals.updateIndex( 1, base+k );
myvals.updateIndex( 2, base+k );
myvals.updateIndex( 3, base+k );
}
}
}
}
Expand Down
Loading