Skip to content

Commit

Permalink
Merge pull request #219 from ji-group/action
Browse files Browse the repository at this point in the history
Multi-treading for action
  • Loading branch information
xji3 authored May 3, 2024
2 parents 2345cf0 + ca3fb7c commit 4d5644e
Show file tree
Hide file tree
Showing 6 changed files with 666 additions and 522 deletions.
1 change: 1 addition & 0 deletions .github/workflows/build.yml
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,7 @@ jobs:
- name: Install (Linux)
if: runner.os == 'Linux' && matrix.name != 'windows'
run: |
sudo apt-get update
sudo apt-get install -y cmake ninja-build ccache pkg-config
sudo apt-get install -y openjdk-17-jdk libeigen3-dev ocl-icd-opencl-dev
Expand Down
29 changes: 24 additions & 5 deletions CMakeLists.txt
Original file line number Diff line number Diff line change
Expand Up @@ -62,7 +62,8 @@ else ()
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -pthread")
endif ()

if (NOT ${CMAKE_BUILD_TYPE} MATCHES Debug)

if (NOT "${CMAKE_BUILD_TYPE}" MATCHES Debug)
if(MSVC)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} /Ox /Ot")
else()
Expand Down Expand Up @@ -109,11 +110,22 @@ if(BUILD_OPENMP)
endif(BUILD_OPENMP)

if(BUILD_ACTION)
find_package (Eigen3 3.3 REQUIRED NO_MODULE)
include(FetchContent)
set(EIGEN_VERSION 3.4)
find_package(Eigen3 ${EIGEN_VERSION} QUIET)
if(NOT EIGEN3_FOUND)
message(" Eigen3: library not found - downloading.")
set(BUILD_TESTING OFF CACHE INTERNAL "")
FetchContent_Declare(eigen
GIT_REPOSITORY https://gitlab.com/libeigen/eigen.git
GIT_TAG ${EIGEN_VERSION}
GIT_SHALLOW ON)
FetchContent_MakeAvailable(eigen)
unset(BUILD_TESTING CACHE)
message(" Eigen3: downloaded.")
endif()

set(CMAKE_CXX_STANDARD 17)
if (NOT TARGET Eigen3::Eigen)
message(FATAL_ERROR "No Eigen3 Found")
endif (NOT TARGET Eigen3::Eigen)
add_definitions(-DBEAGLE_ACTION)
# add_definitions(-DBEAGLE_DEBUG_FLOW)
endif(BUILD_ACTION)
Expand Down Expand Up @@ -180,6 +192,13 @@ if (BEAGLE_OPTIMIZE_FOR_NATIVE_ARCH)
endif()
endif()


include(CheckIncludeFileCXX)
check_include_file_cxx("cpuid.h" HAVE_CPUID_H)
if (HAVE_CPUID_H)
set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -DHAVE_CPUID_H")
endif()

add_subdirectory(libhmsbeagle)

link_libraries(hmsbeagle)
Expand Down
20 changes: 18 additions & 2 deletions examples/actiontest/actiontest.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -125,6 +125,7 @@ void printFlags(long inFlags) {
if (inFlags & BEAGLE_FLAG_VECTOR_SSE) fprintf(stdout, " VECTOR_SSE");
if (inFlags & BEAGLE_FLAG_VECTOR_AVX) fprintf(stdout, " VECTOR_AVX");
if (inFlags & BEAGLE_FLAG_THREADING_NONE) fprintf(stdout, " THREADING_NONE");
if (inFlags & BEAGLE_FLAG_THREADING_CPP) fprintf(stdout, " THREADING_CPP");
if (inFlags & BEAGLE_FLAG_THREADING_OPENMP) fprintf(stdout, " THREADING_OPENMP");
if (inFlags & BEAGLE_FLAG_FRAMEWORK_CPU) fprintf(stdout, " FRAMEWORK_CPU");
if (inFlags & BEAGLE_FLAG_FRAMEWORK_CUDA) fprintf(stdout, " FRAMEWORK_CUDA");
Expand All @@ -147,8 +148,7 @@ int main( int argc, const char* argv[] )
}
fprintf(stdout, "\n");

// bool scaling = true;
bool scaling = false; // disable scaling for now
bool scaling = true;

bool doJC = true;

Expand All @@ -171,6 +171,19 @@ int main( int argc, const char* argv[] )

bool useGpu = argc > 1 && strcmp(argv[1] , "--gpu") == 0;

bool useThreading = false;
for(int i=1;i<argc;i++)
if (!strcmp(argv[i],"--threading"))
useThreading = true;

for(int i=1;i<argc;i++)
if (!strcmp(argv[i],"--help"))
{
std::cerr<<"Flag: --gpu\n";
std::cerr<<"Flag: --threading\n";
std::exit(1);
}

bool useTipStates = false;

int whichDevice = -1;
Expand Down Expand Up @@ -207,6 +220,9 @@ int main( int argc, const char* argv[] )
requirementFlags |= BEAGLE_FLAG_VECTOR_NONE;
}

if (useThreading)
preferenceFlags |= BEAGLE_FLAG_THREADING_CPP;

// create an instance of the BEAGLE library
int instance = beagleCreateInstance(
3, /**< Number of tip data elements (input) */
Expand Down
122 changes: 47 additions & 75 deletions libhmsbeagle/CPU/BeagleCPUActionImpl.h
Original file line number Diff line number Diff line change
Expand Up @@ -82,29 +82,23 @@ namespace beagle {
using BeagleCPUImpl<BEAGLE_CPU_ACTION_DOUBLE>::kFlags;
// SpMatrix** gScaledQs;
int kPartialsCacheOffset;
MapType** gMappedPartials;
MapType** gMappedPartialCache;
// using BeagleCPUImpl<BEAGLE_CPU_ACTION_DOUBLE>::gStateFrequencies;
// using BeagleCPUImpl<BEAGLE_CPU_ACTION_DOUBLE>::gTipStates;
// using BeagleCPUImpl<BEAGLE_CPU_ACTION_DOUBLE>::gTransitionMatrices;
double* gIntegrationTmp;
// double* gLeftPartialTmp;
// double* gRightPartialTmp;
SpMatrix* gInstantaneousMatrices;
SpMatrix* gBs;
double* gMuBs;
double* gB1Norms;
int* gEigenMaps;
double* gEdgeMultipliers;
std::map<int, SpMatrix>* powerMatrices;
std::map<int, double>* ds;
int* gHighestPowers;
std::vector<SpMatrix> gInstantaneousMatrices;
std::vector<SpMatrix> gBs;
std::vector<double> gMuBs;
std::vector<double> gB1Norms;
std::vector<int> gEigenMaps;
std::vector<double> gEdgeMultipliers;
std::vector<std::vector<double>> ds;
SpMatrix identity;
SpMatrix* gScaledQTransposeTmp;
MapType* gMappedIntegrationTmp;
// MapType* gMappedLeftPartialTmp;
// MapType* gMappedRightPartialTmp;
double* gRescaleTmp;
const int mMax = 55;
std::map<int, double> thetaConstants = {
//The first 30 values are from table A.3 of Computing Matrix Functions.
Expand Down Expand Up @@ -169,31 +163,32 @@ namespace beagle {
long long preferenceFlags,
long long requirementFlags);

virtual int setPartials(int bufferIndex,
const double* inPartials);

// virtual int setCategoryRates(const double* inCategoryRates);

virtual int setTipPartials(int tipIndex,
const double* inPartials);
int setTipStates(int tipIndex, const int* inStates);

virtual int updatePartials(const int *operations,
int operationCount,
int cumulativeScalingIndex);
protected:
virtual int upPartials(bool byPartition,
const int *operations,
int operationCount,
int cumulativeScalingIndex);

virtual int updatePrePartials(const int *operations,
int operationCount,
int cumulativeScalingIndex);
// protected:
// virtual int getPaddedPatternsModulus();
virtual int upPrePartials(bool byPartition,
const int *operations,
int operationCount,
int cumulativeScalingIndex);

private:
virtual void rescalePartials(MapType *destP,
double *scaleFactors,
double *cumulativeScaleFactors,
const int fillWithOnes);
inline MapType partialsMap(int index, int category, int startPattern, int endPattern);

inline MapType partialsMap(int index, int category);

inline MapType partialsCacheMap(int index, int category, int startPattern, int endPattern);

inline MapType partialsCacheMap(int index, int category);

// virtual int getPaddedPatternsModulus();

private:
virtual int setEigenDecomposition(int eigenIndex,
const double *inEigenVectors,
const double *inInverseEigenVectors,
Expand All @@ -206,58 +201,35 @@ namespace beagle {
const double* edgeLengths,
int count);

void simpleAction2(MapType *destP, MapType *partials, int edgeIndex,
bool transpose);
void simpleAction2(MapType& destP, MapType& partials, int edgeIndex,
int category, bool transpose) const;

void simpleAction3(MapType& destP, MapType& partials, int edgeIndex,
int category, bool transpose) const;

void simpleAction(MapType* destP,
MapType* partials,
SpMatrix* matrix);
void calcPartialsPartials2(int destPIndex,
int partials1Index,
int edgeIndex1,
int partials2Index,
int edgeIndex2,
int startPattern,
int endPattern);

void calcPartialsPartials(MapType* destP,
MapType* partials1,
SpMatrix* matrices1,
MapType* partials2,
SpMatrix* matrices2);

void calcPartialsPartials2(MapType *destP, MapType *partials1,
MapType *partials2, int edgeIndex1,
int edgeIndex2, MapType *partialCache1,
MapType *partialCache2);

void calcPrePartialsPartials(MapType *destP,
MapType *partials1,
SpMatrix *matrices1,
MapType *partials2,
SpMatrix *matrices2);

void calcPrePartialsPartials2(MapType *destP, MapType *partials1,
MapType *partials2, int edgeIndex1,
int edgeIndex2, MapType *partialCache2);

// Return (m,s)
std::tuple<int,int> getStatistics(double A1Norm,
SpMatrix * matrix,
double t,
int nCol);
void calcPrePartialsPartials2(int destPIndex,
int partials1Index,
int edgeIndex1,
int partials2Index,
int edgeIndex2,
int startPattern,
int endPattern);

// Return (m,s)
std::tuple<int,int> getStatistics2(double t, int nCol, double edgeMultiplier,
int eigenIndex);

double getDValue(int p,
std::map<int, double> &d,
std::map<int, SpMatrix> &powerMatrices);

double getDValue2(int p, int eigenIndex);

double normP1(SpMatrix * matrix);

double normPInf(SpMatrix* matrix);
double normPInf(MapType matrix);
double normPInf(MatrixXd * matrix);
int eigenIndex) const;

double getDValue(int p, int eigenIndex) const;

double getPMax() const;
};


Expand Down
Loading

0 comments on commit 4d5644e

Please sign in to comment.