Merge pull request #20 from saudzahirr/eigen

Add thirdParty libs

CMakeLists.txt

@@ -7,6 +7,15 @@ if(SKBUILD)
   message(STATUS "The project is built using scikit-build")
 endif()
 
+# Find and configure OpenMP
+find_package(OpenMP REQUIRED)
+
+if(OpenMP_CXX_FOUND)
+    message(STATUS "OpenMP found and enabled")
+    add_definitions(-DUSE_OPENMP)
+    set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${OpenMP_CXX_FLAGS}")
+endif()
+
 set(SOURCES
     src/eigenlib/eigen_API.cpp
     src/eigenlib/logger.cpp
@@ -19,13 +28,28 @@ set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${PROJECT_SOURCE_DIR}/src/lib)
 
 include_directories(src/lib)
 
-set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -fPIC -fopenmp")
+# Add subdirectory for third-party libraries
+add_subdirectory(src/eigenlib/thirdParty)
+
+# Include directories for Eigen and fmt
+include_directories(src/eigenlib/thirdParty/Eigen_3.4.0)
+include_directories(src/eigenlib/thirdParty/fmt-10.2.1/include)
 
 configure_file(
   "${PROJECT_SOURCE_DIR}/src/eigenlib/version.h.in"
   "${PROJECT_SOURCE_DIR}/src/eigenlib/version.h"
 )
 
+# Add shared library target for your main project
 add_library(eigen SHARED ${SOURCES})
 
+# Link the fmt static library with your eigen shared library
+target_link_libraries(eigen PRIVATE fmt)
+
+# Link OpenMP if found
+if(OpenMP_CXX_FOUND)
+    target_link_libraries(eigen PRIVATE OpenMP::OpenMP_CXX)
+endif()
+
+# Install the eigen library
 install(TARGETS eigen LIBRARY DESTINATION lib)

pyproject.toml

@@ -5,6 +5,8 @@ requires = [
   "scikit-build",
   "cmake",
   "ninja",
+  "numpy",
+  "scipy",
 ]
 
 [project]
@@ -17,8 +19,8 @@ classifiers = [
   "License :: MIT",
   "Operating System :: Linux",
 ]
-description = "Numerical Algorithm for Eigen Value Problems."
-keywords = ["numerical"]
+description = "Numerical Algorithm for Eigenvalue Problems using Rayleigh Quotient and von Mises Iteration with Deflation."
+keywords = ["numerical", "eigenvalues", "rayleigh", "von_mises"]
 license = {file = "LICENSE"}
 requires-python = ">=3.11"
 

src/eigen.py

@@ -1,6 +1,19 @@
+"""
+Copyright (c) 2024 Saud Zahir
+
+This file is part of vonMises.
+
+vonMises is free software; you can redistribute it and/or
+modify it under the terms of the MIT License.
+
+vonMises is distributed in the hope that it will be useful,
+but WITHOUT ANY WARRANTY; without even the implied warranty of
+MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MIT
+License for more details.
+"""
+
 import ctypes
 import sys
-import time
 from pathlib import Path
 
 if sys.platform == "win32":
@@ -19,20 +32,43 @@ def solve(self, A):
         solveEigenValueProblem = self.lib.solveEigenValueProblem
         solveEigenValueProblem.argtypes = [
             ctypes.POINTER(ctypes.POINTER(ctypes.c_double)),
-            ctypes.c_int,
+            ctypes.c_size_t,
+            ctypes.POINTER(ctypes.c_double),
+            ctypes.POINTER(ctypes.POINTER(ctypes.c_double)),
         ]
         solveEigenValueProblem.restype = None
 
         size = len(A)
 
+        # Allocate and flatten the matrix
         A_ptr = (ctypes.POINTER(ctypes.c_double) * size)()
         for i, row in enumerate(A):
             A_ptr[i] = (ctypes.c_double * len(row))(*row)
 
-        start_time = time.time()
+        # Allocate memory for eigenvalues and eigenvectors
+        eigenvalues = (ctypes.c_double * size)()
+        eigenvectors = (ctypes.POINTER(ctypes.c_double) * size)()
+        for i in range(size):
+            eigenvectors[i] = (ctypes.c_double * size)()
+
+        # Call the C++ function
+        solveEigenValueProblem(A_ptr, size, eigenvalues, eigenvectors)
+
+        eigenvalues = list(eigenvalues)
+        eigenvectors = [list(eigenvectors[i][0:size]) for i in range(size)]
+
+        return eigenvalues, eigenvectors
+
+
+def eigen(A):
+    return EigenValueProblem(eigen_lib).solve(A)
+
+
+def eigenvalues(A):
+    _eigenvalues, _ = eigen(A)
+    return _eigenvalues
 
-        solveEigenValueProblem(A_ptr, size)
 
-        end_time = time.time()
-        execution_time = end_time - start_time
-        print("Execution time:", execution_time, "seconds")
+def eigenvectors(A):
+    _, _eigenvectors = eigen(A)
+    return _eigenvectors

src/eigenlib/eigen_API.cpp

@@ -1,15 +1,48 @@
+/*
+ * Copyright (c) 2024 Saud Zahir
+ *
+ * This file is part of vonMises.
+ *
+ * vonMises is free software; you can redistribute it and/or
+ * modify it under the terms of the MIT License.
+ *
+ * vonMises is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MIT
+ * License for more details.
+ */
+
+
 #include "eigen_API.h"
 
+/**
+ * @file
+ * Eigen value problem solver using deflation method power iteration
+ */
 
-void solveEigenValueProblem(double** A, int size) {
+void solveEigenValueProblem(double** A, size_t size, double* eigenvalues, double** eigenvectors) {
+    INFO("----------------------------------------------");
+    INFO("CMake Version   :: {}", CMake_VERSION);
+    INFO("GCC Version     :: {}", gcc_VERSION);
+    INFO("OMP max threads :: {}", omp_get_max_threads());
+    INFO("----------------------------------------------");
+
+    MatrixXd _A =  convertToEigenMatrix(A, size);
     // Power Iteration with Deflation
-    for (int i = 0; i < size; i++) {
-        double* x = VonMisesIterationMethod(A, size);
-        double c = RayleighQuotient(A, x, size);
-        DEBUG_OUTFILE("*EIGEN_VALUE: " + to_string(c));
-        DEBUG_OUTFILE("*EIGEN_VECTOR: " + getVectorString(x, size));
-        cout << endl;
-        A = matrixSub(A, matrixScalarProduct(outerProduct(x, size, x, size), size, size, c), size, size);
-        delete[] x;
+    for (size_t a = 0; a < size; a++) {
+        double start_time = clocktime();
+        VectorXd x = vonMisesIterationMethod(_A);
+        DEBUG("Execution time for Von Mises Method: {} seconds", deltaTime(start_time));
+        start_time = clocktime();
+        double c = rayleighQuotient(_A, x);
+        DEBUG("Execution time for Rayleigh Quotient Method: {} seconds", deltaTime(start_time));
+        // Store the computed eigenvalue
+        eigenvalues[a] = c;
+
+        // Store the eigenvector
+        for (size_t i = 0; i < size; i++) {
+            eigenvectors[a][i] = x(i);
+        }
+        _A = _A - c * (x * x.transpose());
     }
 }

src/eigenlib/eigen_API.h

@@ -1,3 +1,18 @@
+/*
+ * Copyright (c) 2024 Saud Zahir
+ *
+ * This file is part of vonMises.
+ *
+ * vonMises is free software; you can redistribute it and/or
+ * modify it under the terms of the MIT License.
+ *
+ * vonMises is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the MIT
+ * License for more details.
+ */
+
+
 #ifndef EIGEN_H
 #define EIGEN_H
 
@@ -7,11 +22,14 @@
 #define DLLEXPORT
 #endif
 
+#include <omp.h>
+
 #include "rayleigh_quotient.h"
 #include "von_mises.h"
+#include "version.h"
 
 extern "C" {
-    DLLEXPORT void solveEigenValueProblem(double** A, int size);
+    DLLEXPORT void solveEigenValueProblem(double** A, size_t size, double* eigenvalues, double** eigenvectors);
 }
 
-#endif
+#endif