Multi-treading for action

.github/workflows/build.yml

@@ -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
 

CMakeLists.txt

@@ -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()
@@ -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)
@@ -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)

examples/actiontest/actiontest.cpp

@@ -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");
@@ -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;
 
@@ -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;
@@ -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) */

libhmsbeagle/CPU/BeagleCPUActionImpl.h

@@ -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.
@@ -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,
@@ -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;
         };