Comment cupy functions, add todo list, add tests, add secular function

Author: lorenzotomada

README.md

@@ -3,3 +3,12 @@ Names: Lorenzo Tomada, Gaspare Li Causi
 email: ltomada@sissa.it, glicausi@sissa.it
 
 This repository contains the final project for the course in Development Tools in Scientific Computing.
+
+
+TO DO:
+1) Integrate the new code in the parallel solver
+2) Profile runtime and memory usage, saving the results and plotting
+3) Runtime vs matrix size comparison (follow in detail the instructions on the course repo)
+4) Accuracy vs efficiency
+5) Script to run on Ulysses
+6) Add missing tests

requirements.txt

@@ -1,62 +1,44 @@
 alabaster==1.0.0
-asttokens==3.0.0
 babel==2.17.0
 black==25.1.0
 certifi==2025.1.31
 charset-normalizer==3.4.1
 click==8.1.8
 contourpy==1.3.1
-cupy-cuda12x==13.4.0
 cycler==0.12.1
-decorator==5.2.1
-distro==1.9.0
 docutils==0.21.2
 exceptiongroup==1.2.2
-executing==2.2.0
-fastrlock==0.8.3
 fonttools==4.56.0
 idna==3.10
 imagesize==1.4.1
 iniconfig==2.0.0
-ipython==8.33.0
-ipython_pygments_lexers==1.1.1
-jedi==0.19.2
 Jinja2==3.1.5
 kiwisolver==1.4.8
 line_profiler==4.2.0
 llvmlite==0.44.0
 MarkupSafe==3.0.2
 matplotlib==3.10.0
-matplotlib-inline==0.1.7
 memory-profiler==0.61.0
 mypy-extensions==1.0.0
 numba==0.61.0
 numpy==2.1.3
-nvprof==0.2
 packaging==24.2
 pandas==2.2.3
-parso==0.8.4
 pathspec==0.12.1
-pexpect==4.9.0
 pillow==11.1.0
 platformdirs==4.3.6
 pluggy==1.5.0
-prompt_toolkit==3.0.50
 psutil==7.0.0
-ptyprocess==0.7.0
-pure_eval==0.2.3
 pybind11==2.13.6
 Pygments==2.19.1
 pyparsing==3.2.1
 pytest==8.3.4
 python-dateutil==2.9.0.post0
-pytz==2025.1
+pytz==2025.2
 PyYAML==6.0.2
 requests==2.32.3
-roman-numerals-py==3.1.0
-scikit-build==0.18.1
 scipy==1.15.2
-setuptools==75.8.0
+setuptools==78.1.1
 six==1.17.0
 snowballstemmer==2.2.0
 Sphinx==8.1.3
@@ -68,12 +50,8 @@ sphinxcontrib-jquery==4.1
 sphinxcontrib-jsmath==1.0.1
 sphinxcontrib-qthelp==2.0.0
 sphinxcontrib-serializinghtml==2.0.0
-stack-data==0.6.3
-tokenize_rt==6.1.0
 tomli==2.2.1
-traitlets==5.14.3
 typing_extensions==4.12.2
-tzdata==2025.1
+tzdata==2025.2
 urllib3==2.3.0
-wcwidth==0.2.13
 wheel==0.45.1

scripts/profiling.py

@@ -3,18 +3,18 @@
 from pyclassify import (
     eigenvalues_np,
     eigenvalues_sp,
-    #eigenvalues_cp,
+    # eigenvalues_cp,
     power_method,
     power_method_numba,
-    #power_method_cp,
+    # power_method_cp,
     QR,
-    #QR_cp,
+    # QR_cp,
 )
 from pyclassify.utils import (
     make_symmetric,
     read_config,
     profile_with_cprofile,
-    #profile_with_cupy_profiler,
+    # profile_with_cupy_profiler,
 )
 import numpy as np
 import scipy.sparse as sp
@@ -25,7 +25,7 @@
 
 seed = 8422
 random.seed(seed)
-#cp.random.seed(seed)
+# cp.random.seed(seed)
 np.random.seed(seed)
 
 
@@ -51,7 +51,7 @@
 A = U @ A @ U.T
 A = make_symmetric(A)
 A = sp.csr_matrix(A)
-#A_cp = cpsp.csr_matrix(A)
+# A_cp = cpsp.csr_matrix(A)
 
 
 log_file = "./logs/timings.csv"
@@ -73,6 +73,6 @@
 )
 
 
-#profile_with_cupy_profiler(log_file, dim, "eigenvalues_cp", eigenvalues_cp, A_cp)
-#profile_with_cupy_profiler(log_file, dim, "power_method_cp", power_method_cp, A_cp)
-#profile_with_cupy_profiler(log_file, dim, "QR_cp", QR_cp, A_cp, q0=cp.random.rand(dim), tol=1e-3, max_iter=iteration_factor * dim)
+# profile_with_cupy_profiler(log_file, dim, "eigenvalues_cp", eigenvalues_cp, A_cp)
+# profile_with_cupy_profiler(log_file, dim, "power_method_cp", power_method_cp, A_cp)
+# profile_with_cupy_profiler(log_file, dim, "QR_cp", QR_cp, A_cp, q0=cp.random.rand(dim), tol=1e-3, max_iter=iteration_factor * dim)

scripts/run.py

@@ -1,20 +1,20 @@
-#import cupyx.scipy.sparse as cpsp
-#import cupy as cp
+# import cupyx.scipy.sparse as cpsp
+# import cupy as cp
 from pyclassify import (
     eigenvalues_np,
     eigenvalues_sp,
-    #eigenvalues_cp,
+    # eigenvalues_cp,
     power_method,
     power_method_numba,
-    #power_method_cp,
+    # power_method_cp,
     QR,
-    #QR_cp,
+    # QR_cp,
 )
 from pyclassify.utils import (
     make_symmetric,
     read_config,
     profile_with_cprofile,
-    #profile_with_cupy_profiler,
+    # profile_with_cupy_profiler,
 )
 import numpy as np
 import scipy.sparse as sp
@@ -23,13 +23,13 @@
 import argparse
 
 
-#cp.cuda.Device(0).use()
-#cp.get_default_memory_pool().free_all_blocks()
+# cp.cuda.Device(0).use()
+# cp.get_default_memory_pool().free_all_blocks()
 
 
 seed = 8422
 random.seed(seed)
-#cp.random.seed(seed)
+# cp.random.seed(seed)
 np.random.seed(seed)
 
 
@@ -55,7 +55,7 @@
 A = U @ A @ U.T
 A = make_symmetric(A)
 A = sp.csr_matrix(A)
-#A_cp = cpsp.csr_matrix(A)
+# A_cp = cpsp.csr_matrix(A)
 
 
 log_file = "./logs/timings.csv"
@@ -76,6 +76,6 @@
     log_file, dim, "QR", QR, A.toarray(), max_iter=iteration_factor * dim
 )
 
-#profile_with_cupy_profiler(log_file, dim, "eigenvalues_cp", eigenvalues_cp, A_cp)
-#profile_with_cupy_profiler(log_file, dim, "power_method_cp", power_method_cp, A_cp)
-#profile_with_cupy_profiler(log_file, dim, "QR_cp", QR_cp, A_cp, q0=cp.random.rand(dim), tol=1e-3, max_iter=iteration_factor * dim)
+# profile_with_cupy_profiler(log_file, dim, "eigenvalues_cp", eigenvalues_cp, A_cp)
+# profile_with_cupy_profiler(log_file, dim, "power_method_cp", power_method_cp, A_cp)
+# profile_with_cupy_profiler(log_file, dim, "QR_cp", QR_cp, A_cp, q0=cp.random.rand(dim), tol=1e-3, max_iter=iteration_factor * dim)

setup.py

@@ -41,8 +41,8 @@ def build_extension(self, ext):
 setup(
     name="pyclassify",
     version="0.0.1",
-    author="Lorenzo, Gaspare",
-    author_email="mail.it",
+    author="Lorenzo Tomada, Gaspare Li Causi",
+    author_email="ltomada@sissa.it, glicausi@sissa.it",
     description="Final project",
     long_description="Eigenvalue computation",
     ext_modules=[CMakeExtension("pyclassify.QR_cpp")],

shell/submit.sh

@@ -1,3 +1,4 @@
 #!/bin/bash
 
-pytest -v > output_pytest.txt
+module load cmake/3.29.1
+module load intel/2021.2

src/pyclassify/__init__.py

@@ -1,27 +1,39 @@
 __all__ = [
     "eigenvalues_np",
     "eigenvalues_sp",
-    #"eigenvalues_cp",
+    # "eigenvalues_cp",
     "power_method",
     "power_method_numba",
-    #"power_method_cp",
+    # "power_method_cp",
     "EigenSolver",
-    #"Lanczos_PRO_cp",
-    #"QR_method_cp",
-    #"QR_cp",
+    # "Lanczos_PRO_cp",
+    # "QR_method_cp",
+    # "QR_cp",
 ]
 
 from .QR_cpp import QR_algorithm, Eigen_value_calculator
 
 from .eigenvalues import (
     eigenvalues_np,
     eigenvalues_sp,
-    #eigenvalues_cp,
+    # eigenvalues_cp,
     power_method,
     power_method_numba,
-    #power_method_cp,
+    # power_method_cp,
     EigenSolver,
-    #Lanczos_PRO_cp,
-    #QR_method_cp,
-    #QR_cp,
+    # Lanczos_PRO_cp,
+    # QR_method_cp,
+    # QR_cp,
+)
+
+from .helpers_secular import (
+    inner_outer_eigs,
+    return_secular_f,
+    secular_function,
+    check_is_root,
+    bisection,
+    compute_outer_zero,
+    compute_psi_s,
+    compute_inner_zero,
+    compute_eigenvalues,
 )

src/pyclassify/eigenvalues.py

@@ -2,9 +2,10 @@
 import scipy.sparse as sp
 import scipy.linalg as spla
 from numpy.linalg import eig, eigh
-#import cupy as cp
-#import cupy.linalg as cpla
-#from cupyx.scipy.sparse.linalg import eigsh as eigsh_cp
+
+# import cupy as cp
+# import cupy.linalg as cpla
+# from cupyx.scipy.sparse.linalg import eigsh as eigsh_cp
 from numba import jit, prange
 from .QR_cpp import QR_algorithm, Eigen_value_calculator
 from pyclassify.utils import (
@@ -14,7 +15,6 @@
 )
 
 
-
 def eigenvalues_np(A, symmetric=True):
     """
     Compute the eigenvalues of a square matrix using NumPy's `eig` or `eigh` function.
@@ -43,7 +43,6 @@ def eigenvalues_np(A, symmetric=True):
     return eigenvalues
 
 
-
 def eigenvalues_sp(A, symmetric=True):
     """
     Compute the eigenvalues of a sparse matrix using SciPy's `eigsh` or `eigs` function.
@@ -73,8 +72,7 @@ def eigenvalues_sp(A, symmetric=True):
     return eigenvalues
 
 
-
-#def eigenvalues_cp(A):
+# def eigenvalues_cp(A):
 #    """
 #    Compute the eigenvalues of a sparse matrix using CuPy's `eigsh` function.
 #
@@ -100,7 +98,6 @@ def eigenvalues_sp(A, symmetric=True):
 #    return eigenvalues
 
 
-
 def power_method(A, max_iter=500, tol=1e-4, x=None):
     """
     Compute the dominant eigenvalue of a square matrix using the power method.
@@ -139,7 +136,6 @@ def power_method(A, max_iter=500, tol=1e-4, x=None):
     return x @ A @ x
 
 
-
 @jit(nopython=True, nogil=True, parallel=True)
 def power_method_numba(A, max_iter=500, tol=1e-4, x=None):
     """
@@ -166,7 +162,9 @@ def power_method_numba(A, max_iter=500, tol=1e-4, x=None):
                     because of numba technicalities.
     """
     if A.shape[0] != A.shape[1]:
-        raise ValueError("Matrix must be square!")  # not possible to use the helper function due to the fact that we are using JIT-compilation
+        raise ValueError(
+            "Matrix must be square!"
+        )  # not possible to use the helper function due to the fact that we are using JIT-compilation
     if x is None:
         x = np.random.rand(A.shape[0])
     x /= np.linalg.norm(x)
@@ -186,7 +184,6 @@ def power_method_numba(A, max_iter=500, tol=1e-4, x=None):
     return x @ A @ x
 
 
-
 class EigenSolver:
     """
     Class solving the eigenvalue problem for a given symmetric matrix.
@@ -384,8 +381,7 @@ def eig(self, diag=None, off_diag=None):
         return np.array(eig), Q_triangular @ self.Q.T
 
 
-
-#def power_method_cp(A, max_iter=500, tol=1e-4, x=None):
+# def power_method_cp(A, max_iter=500, tol=1e-4, x=None):
 #    """
 #    Compute the dominant eigenvalue of a square matrix using the power method.
 #    Implemented using cupy.
@@ -425,7 +421,7 @@ def eig(self, diag=None, off_diag=None):
 #    return x @ A @ x
 #
 #
-#def Lanczos_PRO_cp(A, q=None, m=None, tol=1e-8):
+# def Lanczos_PRO_cp(A, q=None, m=None, tol=1e-8):
 #    r"""
 #    Perform the Lanczos algorithm for symmetric matrices.
 #
@@ -492,7 +488,7 @@ def eig(self, diag=None, off_diag=None):
 #    return Q, cp.array(alpha), cp.array(beta[:-1])
 #
 #
-#def QR_method_cp(diag, off_diag, tol=1e-8, max_iter=1000):
+# def QR_method_cp(diag, off_diag, tol=1e-8, max_iter=1000):
 #    """
 #    Compute the eigenvalues of a tridiagonal matrix using the QR algorithm.
 #
@@ -620,7 +616,7 @@ def eig(self, diag=None, off_diag=None):
 #    return diag, Q
 #
 #
-#def QR_cp(A, q0=None, tol=1e-8, max_iter=1000):
+# def QR_cp(A, q0=None, tol=1e-8, max_iter=1000):
 #    """
 #    Compute the eigenvalues of a square matrix using the QR algorithm.
 #    Done using the Lanczos algorithm to compute the tridiagonal matrix and then the QR