Skip to content

GitLab

Commit 787869c3 authored by peastman's avatar peastman
Browse files

Parallelized extraction of inverse matrix for CCMA

parent 848e1790
Hide whitespace changes
Inline Side-by-side
Showing with 145 additions and 306 deletions
platforms/cuda/src/CudaIntegrationUtilities.cpp
View file @ 787869c3
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2009-2014 Stanford University and the Authors. * * Portions copyright (c) 2009-2015 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
#include "openmm/VirtualSite.h" #include "openmm/VirtualSite.h"
#include "quern.h" #include "quern.h"
#include "CudaExpressionUtilities.h" #include "CudaExpressionUtilities.h"
#include "ReferenceCCMAAlgorithm.h"
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
#include <cstdlib> #include <cstdlib>
...@@ -304,157 +305,54 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S ...@@ -304,157 +305,54 @@ CudaIntegrationUtilities::CudaIntegrationUtilities(CudaContext& context, const S
int numCCMA = (int) ccmaConstraints.size(); int numCCMA = (int) ccmaConstraints.size();
if (numCCMA > 0) { if (numCCMA > 0) {
vector<vector<int> > atomConstraints(context.getNumAtoms()); // Record information needed by ReferenceCCMAAlgorithm.
vector<pair<int, int> > refIndices(numCCMA);
vector<RealOpenMM> refDistance(numCCMA);
for (int i = 0; i < numCCMA; i++) { for (int i = 0; i < numCCMA; i++) {
atomConstraints[atom1[ccmaConstraints[i]]].push_back(i); int index = ccmaConstraints[i];
atomConstraints[atom2[ccmaConstraints[i]]].push_back(i); refIndices[i] = make_pair(atom1[index], atom2[index]);
refDistance[i] = distance[index];
} }
vector<vector<int> > linkedConstraints(numCCMA); vector<RealOpenMM> refMasses(numAtoms);
for (unsigned atom = 0; atom < atomConstraints.size(); atom++) { for (int i = 0; i < numAtoms; ++i)
for (unsigned i = 0; i < atomConstraints[atom].size(); i++) refMasses[i] = (RealOpenMM) system.getParticleMass(i);
for (unsigned j = 0; j < i; j++) {
int c1 = atomConstraints[atom][i]; // Look up angles for CCMA.
int c2 = atomConstraints[atom][j];
linkedConstraints[c1].push_back(c2); vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
linkedConstraints[c2].push_back(c1); for (int i = 0; i < system.getNumForces(); i++) {
const HarmonicAngleForce* force = dynamic_cast<const HarmonicAngleForce*>(&system.getForce(i));
if (force != NULL) {
for (int j = 0; j < force->getNumAngles(); j++) {
int atom1, atom2, atom3;
double angle, k;
force->getAngleParameters(j, atom1, atom2, atom3, angle, k);
angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM) angle));
} }
}
int maxLinks = 0;
for (unsigned i = 0; i < linkedConstraints.size(); i++)
maxLinks = max(maxLinks, (int) linkedConstraints[i].size());
int maxAtomConstraints = 0;
for (unsigned i = 0; i < atomConstraints.size(); i++)
maxAtomConstraints = max(maxAtomConstraints, (int) atomConstraints[i].size());
// Compute the constraint coupling matrix
vector<vector<int> > atomAngles(numAtoms);
HarmonicAngleForce const* angleForce = NULL;
for (int i = 0; i < system.getNumForces() && angleForce == NULL; i++)
angleForce = dynamic_cast<HarmonicAngleForce const*>(&system.getForce(i));
if (angleForce != NULL)
for (int i = 0; i < angleForce->getNumAngles(); i++) {
int particle1, particle2, particle3;
double angle, k;
angleForce->getAngleParameters(i, particle1, particle2, particle3, angle, k);
atomAngles[particle2].push_back(i);
}
vector<vector<pair<int, double> > > matrix(numCCMA);
for (int j = 0; j < numCCMA; j++) {
for (int k = 0; k < numCCMA; k++) {
if (j == k) {
matrix[j].push_back(pair<int, double>(j, 1.0));
continue;
}
double scale;
int cj = ccmaConstraints[j];
int ck = ccmaConstraints[k];
int atomj0 = atom1[cj];
int atomj1 = atom2[cj];
int atomk0 = atom1[ck];
int atomk1 = atom2[ck];
int atoma, atomb, atomc;
double imj0 = 1.0/system.getParticleMass(atomj0);
double imj1 = 1.0/system.getParticleMass(atomj1);
if (atomj0 == atomk0) {
atoma = atomj1;
atomb = atomj0;
atomc = atomk1;
scale = imj0/(imj0+imj1);
}
else if (atomj1 == atomk1) {
atoma = atomj0;
atomb = atomj1;
atomc = atomk0;
scale = imj1/(imj0+imj1);
}
else if (atomj0 == atomk1) {
atoma = atomj1;
atomb = atomj0;
atomc = atomk0;
scale = imj0/(imj0+imj1);
}
else if (atomj1 == atomk0) {
atoma = atomj0;
atomb = atomj1;
atomc = atomk1;
scale = imj1/(imj0+imj1);
}
else
continue; // These constraints are not connected.
// Look for a third constraint forming a triangle with these two.
bool foundConstraint = false;
for (int m = 0; m < numCCMA; m++) {
int other = ccmaConstraints[m];
if ((atom1[other] == atoma && atom2[other] == atomc) || (atom1[other] == atomc && atom2[other] == atoma)) {
double d1 = distance[cj];
double d2 = distance[ck];
double d3 = distance[other];
matrix[j].push_back(pair<int, double>(k, scale*(d1*d1+d2*d2-d3*d3)/(2.0*d1*d2)));
foundConstraint = true;
break;
}
}
if (!foundConstraint && angleForce != NULL) {
// We didn't find one, so look for an angle force field term.
const vector<int>& angleCandidates = atomAngles[atomb];
for (vector<int>::const_iterator iter = angleCandidates.begin(); iter != angleCandidates.end(); iter++) {
int particle1, particle2, particle3;
double angle, ka;
angleForce->getAngleParameters(*iter, particle1, particle2, particle3, angle, ka);
if ((particle1 == atoma && particle3 == atomc) || (particle3 == atoma && particle1 == atomc)) {
matrix[j].push_back(pair<int, double>(k, scale*cos(angle)));
break;
}
}
}
}
}
// Invert it using QR.
vector<int> matrixRowStart;
vector<int> matrixColIndex;
vector<double> matrixValue;
for (int i = 0; i < numCCMA; i++) {
matrixRowStart.push_back(matrixValue.size());
for (int j = 0; j < (int) matrix[i].size(); j++) {
pair<int, double> element = matrix[i][j];
matrixColIndex.push_back(element.first);
matrixValue.push_back(element.second);
} }
} }
matrixRowStart.push_back(matrixValue.size());
int *qRowStart, *qColIndex, *rRowStart, *rColIndex; // Create a ReferenceCCMAAlgorithm. It will build and invert the constraint matrix for us.
double *qValue, *rValue;
int result = QUERN_compute_qr(numCCMA, numCCMA, &matrixRowStart[0], &matrixColIndex[0], &matrixValue[0], NULL, ReferenceCCMAAlgorithm ccma(numAtoms, numCCMA, refIndices, refDistance, refMasses, angles, 0.1);
&qRowStart, &qColIndex, &qValue, &rRowStart, &rColIndex, &rValue); vector<vector<pair<int, double> > > matrix = ccma.getMatrix();
vector<double> rhs(numCCMA);
matrix.clear();
matrix.resize(numCCMA);
for (int i = 0; i < numCCMA; i++) {
// Extract column i of the inverse matrix.
for (int j = 0; j < numCCMA; j++)
rhs[j] = (i == j ? 1.0 : 0.0);
result = QUERN_multiply_with_q_transpose(numCCMA, qRowStart, qColIndex, qValue, &rhs[0]);
result = QUERN_solve_with_r(numCCMA, rRowStart, rColIndex, rValue, &rhs[0], &rhs[0]);
for (int j = 0; j < numCCMA; j++) {
double value = rhs[j]*distance[ccmaConstraints[i]]/distance[ccmaConstraints[j]];
if (abs(value) > 0.1)
matrix[j].push_back(pair<int, double>(i, value));
}
}
QUERN_free_result(qRowStart, qColIndex, qValue);
QUERN_free_result(rRowStart, rColIndex, rValue);
int maxRowElements = 0; int maxRowElements = 0;
for (unsigned i = 0; i < matrix.size(); i++) for (unsigned i = 0; i < matrix.size(); i++)
maxRowElements = max(maxRowElements, (int) matrix[i].size()); maxRowElements = max(maxRowElements, (int) matrix[i].size());
maxRowElements++; maxRowElements++;
// Build the list of constraints for each atom.
vector<vector<int> > atomConstraints(context.getNumAtoms());
for (int i = 0; i < numCCMA; i++) {
atomConstraints[atom1[ccmaConstraints[i]]].push_back(i);
atomConstraints[atom2[ccmaConstraints[i]]].push_back(i);
}
int maxAtomConstraints = 0;
for (unsigned i = 0; i < atomConstraints.size(); i++)
maxAtomConstraints = max(maxAtomConstraints, (int) atomConstraints[i].size());
// Sort the constraints. // Sort the constraints.
vector<int> constraintOrder(numCCMA); vector<int> constraintOrder(numCCMA);
......
platforms/opencl/src/OpenCLIntegrationUtilities.cpp
View file @ 787869c3
...@@ -6,7 +6,7 @@ ...@@ -6,7 +6,7 @@
* Biological Structures at Stanford, funded under the NIH Roadmap for * * Biological Structures at Stanford, funded under the NIH Roadmap for *
* Medical Research, grant U54 GM072970. See https://simtk.org. * * Medical Research, grant U54 GM072970. See https://simtk.org. *
* * * *
* Portions copyright (c) 2009-2014 Stanford University and the Authors. * * Portions copyright (c) 2009-2015 Stanford University and the Authors. *
* Authors: Peter Eastman * * Authors: Peter Eastman *
* Contributors: * * Contributors: *
* * * *
...@@ -32,6 +32,7 @@ ...@@ -32,6 +32,7 @@
#include "openmm/VirtualSite.h" #include "openmm/VirtualSite.h"
#include "quern.h" #include "quern.h"
#include "OpenCLExpressionUtilities.h" #include "OpenCLExpressionUtilities.h"
#include "ReferenceCCMAAlgorithm.h"
#include <algorithm> #include <algorithm>
#include <cmath> #include <cmath>
#include <cstdlib> #include <cstdlib>
...@@ -323,157 +324,54 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c ...@@ -323,157 +324,54 @@ OpenCLIntegrationUtilities::OpenCLIntegrationUtilities(OpenCLContext& context, c
int numCCMA = (int) ccmaConstraints.size(); int numCCMA = (int) ccmaConstraints.size();
if (numCCMA > 0) { if (numCCMA > 0) {
vector<vector<int> > atomConstraints(context.getNumAtoms()); // Record information needed by ReferenceCCMAAlgorithm.
vector<pair<int, int> > refIndices(numCCMA);
vector<RealOpenMM> refDistance(numCCMA);
for (int i = 0; i < numCCMA; i++) { for (int i = 0; i < numCCMA; i++) {
atomConstraints[atom1[ccmaConstraints[i]]].push_back(i); int index = ccmaConstraints[i];
atomConstraints[atom2[ccmaConstraints[i]]].push_back(i); refIndices[i] = make_pair(atom1[index], atom2[index]);
} refDistance[i] = distance[index];
vector<vector<int> > linkedConstraints(numCCMA);
for (unsigned atom = 0; atom < atomConstraints.size(); atom++) {
for (unsigned i = 0; i < atomConstraints[atom].size(); i++)
for (unsigned j = 0; j < i; j++) {
int c1 = atomConstraints[atom][i];
int c2 = atomConstraints[atom][j];
linkedConstraints[c1].push_back(c2);
linkedConstraints[c2].push_back(c1);
}
} }
int maxLinks = 0; vector<RealOpenMM> refMasses(numAtoms);
for (unsigned i = 0; i < linkedConstraints.size(); i++) for (int i = 0; i < numAtoms; ++i)
maxLinks = max(maxLinks, (int) linkedConstraints[i].size()); refMasses[i] = (RealOpenMM) system.getParticleMass(i);
int maxAtomConstraints = 0;
for (unsigned i = 0; i < atomConstraints.size(); i++)
maxAtomConstraints = max(maxAtomConstraints, (int) atomConstraints[i].size());
// Compute the constraint coupling matrix
vector<vector<int> > atomAngles(numAtoms);
HarmonicAngleForce const* angleForce = NULL;
for (int i = 0; i < system.getNumForces() && angleForce == NULL; i++)
angleForce = dynamic_cast<HarmonicAngleForce const*>(&system.getForce(i));
if (angleForce != NULL)
for (int i = 0; i < angleForce->getNumAngles(); i++) {
int particle1, particle2, particle3;
double angle, k;
angleForce->getAngleParameters(i, particle1, particle2, particle3, angle, k);
atomAngles[particle2].push_back(i);
}
vector<vector<pair<int, double> > > matrix(numCCMA);
for (int j = 0; j < numCCMA; j++) {
for (int k = 0; k < numCCMA; k++) {
if (j == k) {
matrix[j].push_back(pair<int, double>(j, 1.0));
continue;
}
double scale;
int cj = ccmaConstraints[j];
int ck = ccmaConstraints[k];
int atomj0 = atom1[cj];
int atomj1 = atom2[cj];
int atomk0 = atom1[ck];
int atomk1 = atom2[ck];
int atoma, atomb, atomc;
double imj0 = 1.0/system.getParticleMass(atomj0);
double imj1 = 1.0/system.getParticleMass(atomj1);
if (atomj0 == atomk0) {
atoma = atomj1;
atomb = atomj0;
atomc = atomk1;
scale = imj0/(imj0+imj1);
}
else if (atomj1 == atomk1) {
atoma = atomj0;
atomb = atomj1;
atomc = atomk0;
scale = imj1/(imj0+imj1);
}
else if (atomj0 == atomk1) {
atoma = atomj1;
atomb = atomj0;
atomc = atomk0;
scale = imj0/(imj0+imj1);
}
else if (atomj1 == atomk0) {
atoma = atomj0;
atomb = atomj1;
atomc = atomk1;
scale = imj1/(imj0+imj1);
}
else
continue; // These constraints are not connected.
// Look for a third constraint forming a triangle with these two.
bool foundConstraint = false; // Look up angles for CCMA.
for (int m = 0; m < numCCMA; m++) {
int other = ccmaConstraints[m]; vector<ReferenceCCMAAlgorithm::AngleInfo> angles;
if ((atom1[other] == atoma && atom2[other] == atomc) || (atom1[other] == atomc && atom2[other] == atoma)) { for (int i = 0; i < system.getNumForces(); i++) {
double d1 = distance[cj]; const HarmonicAngleForce* force = dynamic_cast<const HarmonicAngleForce*>(&system.getForce(i));
double d2 = distance[ck]; if (force != NULL) {
double d3 = distance[other]; for (int j = 0; j < force->getNumAngles(); j++) {
matrix[j].push_back(pair<int, double>(k, scale*(d1*d1+d2*d2-d3*d3)/(2.0*d1*d2))); int atom1, atom2, atom3;
foundConstraint = true; double angle, k;
break; force->getAngleParameters(j, atom1, atom2, atom3, angle, k);
} angles.push_back(ReferenceCCMAAlgorithm::AngleInfo(atom1, atom2, atom3, (RealOpenMM) angle));
}
if (!foundConstraint && angleForce != NULL) {
// We didn't find one, so look for an angle force field term.
const vector<int>& angleCandidates = atomAngles[atomb];
for (vector<int>::const_iterator iter = angleCandidates.begin(); iter != angleCandidates.end(); iter++) {
int particle1, particle2, particle3;
double angle, ka;
angleForce->getAngleParameters(*iter, particle1, particle2, particle3, angle, ka);
if ((particle1 == atoma && particle3 == atomc) || (particle3 == atoma && particle1 == atomc)) {
matrix[j].push_back(pair<int, double>(k, scale*cos(angle)));
break;
}
}
} }
} }
} }
// Invert it using QR. // Create a ReferenceCCMAAlgorithm. It will build and invert the constraint matrix for us.
vector<int> matrixRowStart; ReferenceCCMAAlgorithm ccma(numAtoms, numCCMA, refIndices, refDistance, refMasses, angles, 0.1);
vector<int> matrixColIndex; vector<vector<pair<int, double> > > matrix = ccma.getMatrix();
vector<double> matrixValue;
for (int i = 0; i < numCCMA; i++) {
matrixRowStart.push_back(matrixValue.size());
for (int j = 0; j < (int) matrix[i].size(); j++) {
pair<int, double> element = matrix[i][j];
matrixColIndex.push_back(element.first);
matrixValue.push_back(element.second);
}
}
matrixRowStart.push_back(matrixValue.size());
int *qRowStart, *qColIndex, *rRowStart, *rColIndex;
double *qValue, *rValue;
int result = QUERN_compute_qr(numCCMA, numCCMA, &matrixRowStart[0], &matrixColIndex[0], &matrixValue[0], NULL,
&qRowStart, &qColIndex, &qValue, &rRowStart, &rColIndex, &rValue);
vector<double> rhs(numCCMA);
matrix.clear();
matrix.resize(numCCMA);
for (int i = 0; i < numCCMA; i++) {
// Extract column i of the inverse matrix.
for (int j = 0; j < numCCMA; j++)
rhs[j] = (i == j ? 1.0 : 0.0);
result = QUERN_multiply_with_q_transpose(numCCMA, qRowStart, qColIndex, qValue, &rhs[0]);
result = QUERN_solve_with_r(numCCMA, rRowStart, rColIndex, rValue, &rhs[0], &rhs[0]);
for (int j = 0; j < numCCMA; j++) {
double value = rhs[j]*distance[ccmaConstraints[i]]/distance[ccmaConstraints[j]];
if (abs(value) > 0.1)
matrix[j].push_back(pair<int, double>(i, value));
}
}
QUERN_free_result(qRowStart, qColIndex, qValue);
QUERN_free_result(rRowStart, rColIndex, rValue);
int maxRowElements = 0; int maxRowElements = 0;
for (unsigned i = 0; i < matrix.size(); i++) for (unsigned i = 0; i < matrix.size(); i++)
maxRowElements = max(maxRowElements, (int) matrix[i].size()); maxRowElements = max(maxRowElements, (int) matrix[i].size());
maxRowElements++; maxRowElements++;
// Build the list of constraints for each atom.
vector<vector<int> > atomConstraints(context.getNumAtoms());
for (int i = 0; i < numCCMA; i++) {
atomConstraints[atom1[ccmaConstraints[i]]].push_back(i);
atomConstraints[atom2[ccmaConstraints[i]]].push_back(i);
}
int maxAtomConstraints = 0;
for (unsigned i = 0; i < atomConstraints.size(); i++)
maxAtomConstraints = max(maxAtomConstraints, (int) atomConstraints[i].size());
// Sort the constraints. // Sort the constraints.
vector<int> constraintOrder(numCCMA); vector<int> constraintOrder(numCCMA);
......
platforms/opencl/tests/TestOpenCLVirtualSites.cpp
View file @ 787869c3
...@@ -413,7 +413,7 @@ void testConservationLaws() { ...@@ -413,7 +413,7 @@ void testConservationLaws() {
if (i == 0) if (i == 0)
initialAngularMomentum = angularMomentum; initialAngularMomentum = angularMomentum;
else else
ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.03); ASSERT_EQUAL_VEC(initialAngularMomentum, angularMomentum, 0.05);
integrator.step(1); integrator.step(1);
} }
} }
......
platforms/reference/include/ReferenceCCMAAlgorithm.h
View file @ 787869c3
/* Portions copyright (c) 2006-2013 Stanford University and Simbios. /* Portions copyright (c) 2006-2015 Stanford University and Simbios.
* Contributors: Pande Group * Contributors: Pande Group
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -37,6 +37,7 @@ class OPENMM_EXPORT ReferenceCCMAAlgorithm : public ReferenceConstraintAlgorithm ...@@ -37,6 +37,7 @@ class OPENMM_EXPORT ReferenceCCMAAlgorithm : public ReferenceConstraintAlgorithm
protected: protected:
int _maximumNumberOfIterations; int _maximumNumberOfIterations;
RealOpenMM _elementCutoff;
int _numberOfConstraints; int _numberOfConstraints;
std::vector<std::pair<int, int> > _atomIndices; std::vector<std::pair<int, int> > _atomIndices;
...@@ -66,8 +67,9 @@ public: ...@@ -66,8 +67,9 @@ public:
* @param distance distances for constraints * @param distance distances for constraints
* @param masses atom masses * @param masses atom masses
* @param angles angle force field terms * @param angles angle force field terms
* @param elementCutoff the cutoff for which elements of the inverse matrix to keep
*/ */
ReferenceCCMAAlgorithm(int numberOfAtoms, int numberOfConstraints, const std::vector<std::pair<int, int> >& atomIndices, const std::vector<RealOpenMM>& distance, std::vector<RealOpenMM>& masses, std::vector<AngleInfo>& angles); ReferenceCCMAAlgorithm(int numberOfAtoms, int numberOfConstraints, const std::vector<std::pair<int, int> >& atomIndices, const std::vector<RealOpenMM>& distance, std::vector<RealOpenMM>& masses, std::vector<AngleInfo>& angles, RealOpenMM elementCutoff);
~ReferenceCCMAAlgorithm(); ~ReferenceCCMAAlgorithm();
...@@ -107,6 +109,13 @@ public: ...@@ -107,6 +109,13 @@ public:
*/ */
void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates, void applyToVelocities(std::vector<OpenMM::RealVec>& atomCoordinates,
std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance); std::vector<OpenMM::RealVec>& velocities, std::vector<RealOpenMM>& inverseMasses, RealOpenMM tolerance);
/**
* Get the inverse constraint matrix. Each element represents one column, and contains a list
* of all non-zero elements in the form (index, value).
*/
const std::vector<std::vector<std::pair<int, RealOpenMM> > >& getMatrix() const;
}; };
class ReferenceCCMAAlgorithm::AngleInfo class ReferenceCCMAAlgorithm::AngleInfo
......
platforms/reference/src/SimTKReference/ReferenceCCMAAlgorithm.cpp
View file @ 787869c3
/* Portions copyright (c) 2006-2013 Stanford University and Simbios. /* Portions copyright (c) 2006-2015 Stanford University and Simbios.
* Contributors: Peter Eastman, Pande Group * Contributors: Peter Eastman, Pande Group
* *
* Permission is hereby granted, free of charge, to any person obtaining * Permission is hereby granted, free of charge, to any person obtaining
...@@ -31,6 +31,7 @@ ...@@ -31,6 +31,7 @@
#include "quern.h" #include "quern.h"
#include "openmm/OpenMMException.h" #include "openmm/OpenMMException.h"
#include "openmm/Vec3.h" #include "openmm/Vec3.h"
#include "openmm/internal/ThreadPool.h"
#include <map> #include <map>
using std::map; using std::map;
...@@ -39,13 +40,51 @@ using std::vector; ...@@ -39,13 +40,51 @@ using std::vector;
using std::set; using std::set;
using namespace OpenMM; using namespace OpenMM;
// This class extracts columns from the inverse matrix one at a time. It is done in parallel,
// since this can be very slow.
class ExtractMatrixTask : public ThreadPool::Task {
public:
ExtractMatrixTask(int numConstraints, vector<vector<pair<int, RealOpenMM> > >& matrix, const vector<RealOpenMM>& distance, RealOpenMM elementCutoff,
const int* qRowStart, const int* qColIndex, const int* rRowStart, const int* rColIndex, const double* qValue, const double* rValue) :
numConstraints(numConstraints), matrix(matrix), distance(distance), elementCutoff(elementCutoff), qRowStart(qRowStart), qColIndex(qColIndex),
rRowStart(rRowStart), rColIndex(rColIndex), qValue(qValue), rValue(rValue) {
}
void execute(ThreadPool& pool, int threadIndex) {
vector<double> rhs(numConstraints);
for (int i = threadIndex; i < numConstraints; i += pool.getNumThreads()) {
// Extract column i of the inverse matrix.
for (int j = 0; j < numConstraints; j++)
rhs[j] = (i == j ? 1.0 : 0.0);
QUERN_multiply_with_q_transpose(numConstraints, qRowStart, qColIndex, qValue, &rhs[0]);
QUERN_solve_with_r(numConstraints, rRowStart, rColIndex, rValue, &rhs[0], &rhs[0]);
for (int j = 0; j < numConstraints; j++) {
double value = rhs[j]*distance[j]/distance[i];
if (FABS((RealOpenMM) value) > elementCutoff)
matrix[i].push_back(pair<int, RealOpenMM>(j, (RealOpenMM) value));
}
}
}
private:
int numConstraints;
vector<vector<pair<int, RealOpenMM> > >& matrix;
const vector<RealOpenMM>& distance;
RealOpenMM elementCutoff;
const int *qRowStart, *qColIndex, *rRowStart, *rColIndex;
const double *qValue, *rValue;
};
ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms, ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms,
int numberOfConstraints, int numberOfConstraints,
const vector<pair<int, int> >& atomIndices, const vector<pair<int, int> >& atomIndices,
const vector<RealOpenMM>& distance, const vector<RealOpenMM>& distance,
vector<RealOpenMM>& masses, vector<RealOpenMM>& masses,
vector<AngleInfo>& angles) { vector<AngleInfo>& angles,
RealOpenMM elementCutoff) {
_numberOfConstraints = numberOfConstraints; _numberOfConstraints = numberOfConstraints;
_elementCutoff = elementCutoff;
_atomIndices = atomIndices; _atomIndices = atomIndices;
_distance = distance; _distance = distance;
...@@ -157,19 +196,10 @@ ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms, ...@@ -157,19 +196,10 @@ ReferenceCCMAAlgorithm::ReferenceCCMAAlgorithm(int numberOfAtoms,
&qRowStart, &qColIndex, &qValue, &rRowStart, &rColIndex, &rValue); &qRowStart, &qColIndex, &qValue, &rRowStart, &rColIndex, &rValue);
vector<double> rhs(numberOfConstraints); vector<double> rhs(numberOfConstraints);
_matrix.resize(numberOfConstraints); _matrix.resize(numberOfConstraints);
for (int i = 0; i < numberOfConstraints; i++) { ThreadPool threads;
// Extract column i of the inverse matrix. ExtractMatrixTask task(numberOfConstraints, _matrix, _distance, _elementCutoff, qRowStart, qColIndex, rRowStart, rColIndex, qValue, rValue);
threads.execute(task);
for (int j = 0; j < numberOfConstraints; j++) threads.waitForThreads();
rhs[j] = (i == j ? 1.0 : 0.0);
QUERN_multiply_with_q_transpose(numberOfConstraints, qRowStart, qColIndex, qValue, &rhs[0]);
QUERN_solve_with_r(numberOfConstraints, rRowStart, rColIndex, rValue, &rhs[0], &rhs[0]);
for (int j = 0; j < numberOfConstraints; j++) {
double value = rhs[j]*_distance[i]/_distance[j];
if (FABS((RealOpenMM)value) > 0.02)
_matrix[j].push_back(pair<int, RealOpenMM>(i, (RealOpenMM) value));
}
}
QUERN_free_result(qRowStart, qColIndex, qValue); QUERN_free_result(qRowStart, qColIndex, qValue);
QUERN_free_result(rRowStart, rColIndex, rValue); QUERN_free_result(rRowStart, rColIndex, rValue);
} }
...@@ -290,3 +320,7 @@ void ReferenceCCMAAlgorithm::applyConstraints(vector<RealVec>& atomCoordinates, ...@@ -290,3 +320,7 @@ void ReferenceCCMAAlgorithm::applyConstraints(vector<RealVec>& atomCoordinates,
} }
} }
} }
const vector<vector<pair<int, RealOpenMM> > >& ReferenceCCMAAlgorithm::getMatrix() const {
return _matrix;
}
platforms/reference/src/SimTKReference/ReferenceConstraints.cpp
View file @ 787869c3
...@@ -180,7 +180,7 @@ ReferenceConstraints::ReferenceConstraints(const System& system) : ccma(NULL), s ...@@ -180,7 +180,7 @@ ReferenceConstraints::ReferenceConstraints(const System& system) : ccma(NULL), s
// Create the CCMA object. // Create the CCMA object.
ccma = new ReferenceCCMAAlgorithm(numParticles, numCCMA, ccmaIndices, ccmaDistance, masses, angles); ccma = new ReferenceCCMAAlgorithm(numParticles, numCCMA, ccmaIndices, ccmaDistance, masses, angles, 0.02);
} }
} }
......
Markdown is supported
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