From 6859fa8fc8b528f101444a5108be2ddba2a4e2d1 Mon Sep 17 00:00:00 2001 From: JKutt Date: Fri, 11 Feb 2022 16:40:11 -0800 Subject: [PATCH 1/7] adding in handle sotre and restore functionality --- pydiso/mkl_solver.pyx | 20 ++++++++++++++++++++ 1 file changed, 20 insertions(+) diff --git a/pydiso/mkl_solver.pyx b/pydiso/mkl_solver.pyx index 82a9179..ccdd8e7 100644 --- a/pydiso/mkl_solver.pyx +++ b/pydiso/mkl_solver.pyx @@ -30,6 +30,7 @@ cdef extern from 'mkl.h': void mkl_get_version(MKLVersion* pv) void mkl_set_num_threads(int nth) + int mkl_domain_set_num_threads(int nt, int domain) int mkl_get_max_threads() int mkl_domain_get_max_threads(int domain) @@ -39,6 +40,10 @@ cdef extern from 'mkl.h': ctypedef void * _MKL_DSS_HANDLE_t + void pardiso_handle_store(_MKL_DSS_HANDLE_t pt, char *dirname, int *err) + + void pardiso_handle_restore(_MKL_DSS_HANDLE_t pt, char *dirname, int *err) + void pardiso(_MKL_DSS_HANDLE_t, const int*, const int*, const int*, const int *, const int *, const void *, const int *, const int *, int *, const int *, int *, @@ -177,6 +182,7 @@ ctypedef fused _par_params: cdef class MKLPardisoSolver: cdef _MKL_DSS_HANDLE_t handle[64] + cdef _MKL_DSS_HANDLE_t _MKL_pt cdef _PardisoParams _par cdef _PardisoParams64 _par64 cdef int_t _is_32 @@ -422,6 +428,20 @@ cdef class MKLPardisoSolver: else: self._par.iparm[i] = val + def pardiso_store(dirname): + """ + turns on storing of the factorization files + """ + + pardiso_handle_store(self._MKL_DSS_HANDLE_t) + + def pardiso_restore(dirname): + """ + reads in the stored factorization file + """ + + pass + @property def nnz(self): return self.iparm[17] From baf2a353a59b17449f337289ee73aad4482bf182 Mon Sep 17 00:00:00 2001 From: Juan Date: Sat, 12 Feb 2022 14:55:18 -0800 Subject: [PATCH 2/7] functionality for saving and using stored factorizations --- pydiso/mkl_solver.pyx | 75 +++++++++++++++++++++++++++++++++---------- 1 file changed, 58 insertions(+), 17 deletions(-) diff --git a/pydiso/mkl_solver.pyx b/pydiso/mkl_solver.pyx index ccdd8e7..89f9bb5 100644 --- a/pydiso/mkl_solver.pyx +++ b/pydiso/mkl_solver.pyx @@ -182,7 +182,6 @@ ctypedef fused _par_params: cdef class MKLPardisoSolver: cdef _MKL_DSS_HANDLE_t handle[64] - cdef _MKL_DSS_HANDLE_t _MKL_pt cdef _PardisoParams _par cdef _PardisoParams64 _par64 cdef int_t _is_32 @@ -311,6 +310,8 @@ cdef class MKLPardisoSolver: self._set_A(A.data) self._analyze() self._factored = False + self._store = False + self._flag_dir = "" if factor: self._factor() @@ -428,19 +429,10 @@ cdef class MKLPardisoSolver: else: self._par.iparm[i] = val - def pardiso_store(dirname): - """ - turns on storing of the factorization files - """ + def store_factorization(self, directory="./"): - pardiso_handle_store(self._MKL_DSS_HANDLE_t) - - def pardiso_restore(dirname): - """ - reads in the stored factorization file - """ - - pass + self._store = True + self._flag_dir = directory @property def nnz(self): @@ -535,11 +527,47 @@ cdef class MKLPardisoSolver: cdef _factor(self): #phase = 22 self._factored = False + + if self._store: + try: - err = self._run_pardiso(22) - if err!=0: - raise PardisoError("Factor step error, "+_err_messages[err]) - self._factored = True + flag_file = self._flag_dir + 'flagfile.txt' + + cdef char* call_flag_dir = self._flag_dir + + with open(flag_file, 'x') as f: + f.write('inversion in progress') + + err = self._run_pardiso(22) + + self._pardiso_store(call_flag_dir) + + done_file = self._flag_dir + 'factorization_done.txt' + + with open(done_file, 'w') as f2: + f2.write('done') + + self._factored = True + return + + except FileExistsError: + + # flag file exists, wait for "done" file and read in factorization + cdef char* done_file = self._flag_dir + 'factorization_done.txt' + + while not os.path.isfile(done_file): + time.sleep(1) + + # now read in the factorization from the file + cdef char* call_flag_dir = self._flag_dir + self._pardiso_restore(call_flag_dir) + + else: + + err = self._run_pardiso(22) + if err!=0: + raise PardisoError("Factor step error, "+_err_messages[err]) + self._factored = True cdef _solve(self, void* b, void* x, int_t nrhs_in): #phase = 33 @@ -564,3 +592,16 @@ cdef class MKLPardisoSolver: &phase64, &self._par64.n, self.a, &self._par64.ia[0], &self._par64.ja[0], &self._par64.perm[0], &nrhs64, self._par64.iparm, &self._par64.msglvl, b, x, &error64) return error64 + + cdef _pardiso_store(char *dir_name): + + cdef int_t error=0 + + pardiso_handle_store(self.handle, dir_name, &error) + + cdef _pardiso_restore(char *dirname): + + cdef int_t error=0 + + pardiso_handle_restore(self.handle, dir_name, &error) + From 59e98ea165f59eb08f108abc700729186b8781a2 Mon Sep 17 00:00:00 2001 From: JKutt Date: Tue, 5 Apr 2022 14:51:47 -0700 Subject: [PATCH 3/7] fixed bad syntax causing compilation errors --- pydiso/mkl_solver.pyx | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/pydiso/mkl_solver.pyx b/pydiso/mkl_solver.pyx index 89f9bb5..bba6601 100644 --- a/pydiso/mkl_solver.pyx +++ b/pydiso/mkl_solver.pyx @@ -4,6 +4,7 @@ cimport numpy as np from cython cimport numeric import warnings +import time.time as time import numpy as np import scipy.sparse as sp import os @@ -189,6 +190,7 @@ cdef class MKLPardisoSolver: cdef int_t _factored cdef size_t shape[2] cdef int_t _initialized + cdef char* call_flag_dir cdef void * a @@ -533,14 +535,14 @@ cdef class MKLPardisoSolver: flag_file = self._flag_dir + 'flagfile.txt' - cdef char* call_flag_dir = self._flag_dir + self.call_flag_dir = self._flag_dir with open(flag_file, 'x') as f: f.write('inversion in progress') err = self._run_pardiso(22) - self._pardiso_store(call_flag_dir) + self._pardiso_store(self.call_flag_dir) done_file = self._flag_dir + 'factorization_done.txt' @@ -553,14 +555,14 @@ cdef class MKLPardisoSolver: except FileExistsError: # flag file exists, wait for "done" file and read in factorization - cdef char* done_file = self._flag_dir + 'factorization_done.txt' + done_file = self._flag_dir + 'factorization_done.txt' while not os.path.isfile(done_file): time.sleep(1) # now read in the factorization from the file - cdef char* call_flag_dir = self._flag_dir - self._pardiso_restore(call_flag_dir) + self.call_flag_dir = self._flag_dir + self._pardiso_restore(self.call_flag_dir) else: @@ -593,13 +595,13 @@ cdef class MKLPardisoSolver: &self._par64.perm[0], &nrhs64, self._par64.iparm, &self._par64.msglvl, b, x, &error64) return error64 - cdef _pardiso_store(char *dir_name): + cdef _pardiso_store(self, char *dir_name): cdef int_t error=0 pardiso_handle_store(self.handle, dir_name, &error) - cdef _pardiso_restore(char *dirname): + cdef _pardiso_restore(self, char *dir_name): cdef int_t error=0 From 3dbfcde36a6a0e7142677f6c63c07acfecae1cbd Mon Sep 17 00:00:00 2001 From: JKutt Date: Tue, 5 Apr 2022 15:06:40 -0700 Subject: [PATCH 4/7] fixed bad import of time module --- pydiso/mkl_solver.pyx | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pydiso/mkl_solver.pyx b/pydiso/mkl_solver.pyx index bba6601..cef6300 100644 --- a/pydiso/mkl_solver.pyx +++ b/pydiso/mkl_solver.pyx @@ -4,7 +4,7 @@ cimport numpy as np from cython cimport numeric import warnings -import time.time as time +from time import time import numpy as np import scipy.sparse as sp import os From 11867a14f0acb99aadc9da5b9bfa4d7cc5682393 Mon Sep 17 00:00:00 2001 From: JKutt Date: Tue, 5 Apr 2022 15:45:56 -0700 Subject: [PATCH 5/7] fixed additional cdef variables --- pydiso/mkl_solver.pyx | 2 ++ 1 file changed, 2 insertions(+) diff --git a/pydiso/mkl_solver.pyx b/pydiso/mkl_solver.pyx index cef6300..80d11f0 100644 --- a/pydiso/mkl_solver.pyx +++ b/pydiso/mkl_solver.pyx @@ -191,6 +191,8 @@ cdef class MKLPardisoSolver: cdef size_t shape[2] cdef int_t _initialized cdef char* call_flag_dir + cdef char* _flag_dir + cdef int_t _store cdef void * a From 65fe61624a15be8145a9a349596050d8462b79a5 Mon Sep 17 00:00:00 2001 From: JKutt Date: Thu, 14 Apr 2022 16:26:51 -0700 Subject: [PATCH 6/7] added test for storing factorization file and cleaned up syntax --- pydiso/mkl_solver.pyx | 33 ++++++++++++++++++++++++++------- tests/test.py | 38 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 64 insertions(+), 7 deletions(-) diff --git a/pydiso/mkl_solver.pyx b/pydiso/mkl_solver.pyx index 80d11f0..751c524 100644 --- a/pydiso/mkl_solver.pyx +++ b/pydiso/mkl_solver.pyx @@ -199,7 +199,7 @@ cdef class MKLPardisoSolver: cdef object _data_type cdef object _Adata #a reference to make sure the pointer "a" doesn't get destroyed - def __init__(self, A, matrix_type=None, factor=True, verbose=False): + def __init__(self, A, matrix_type=None, factor=True, verbose=False, store_factorization_dir=None): '''ParidsoSolver(A, matrix_type=None, factor=True, verbose=False) An interface to the intel MKL pardiso sparse matrix solver. @@ -314,8 +314,27 @@ cdef class MKLPardisoSolver: self._set_A(A.data) self._analyze() self._factored = False - self._store = False - self._flag_dir = "" + + # check if we want to store the factorization + if store_factorization_dir is not None: + + # check if the flag files exist. If so delete them so factorization file get overwritten + check_file = store_factorization_dir + 'factorization_done.txt' + + if os.path.exists(check_file): + + second_file_to_remove = store_factorization_dir + "flagfile.txt" + os.remove(check_file) + os.remove(second_file_to_remove) + + self._store = True + flag_dir_ = bytes(store_factorization_dir, 'utf-8') + self._flag_dir = flag_dir_ + + else: + + self._store = False + if factor: self._factor() @@ -433,7 +452,7 @@ cdef class MKLPardisoSolver: else: self._par.iparm[i] = val - def store_factorization(self, directory="./"): + def store_factorization(self, directory=b'./'): self._store = True self._flag_dir = directory @@ -535,7 +554,7 @@ cdef class MKLPardisoSolver: if self._store: try: - flag_file = self._flag_dir + 'flagfile.txt' + flag_file = self._flag_dir.decode("utf-8") + 'flagfile.txt' self.call_flag_dir = self._flag_dir @@ -546,7 +565,7 @@ cdef class MKLPardisoSolver: self._pardiso_store(self.call_flag_dir) - done_file = self._flag_dir + 'factorization_done.txt' + done_file = self._flag_dir.decode("utf-8") + 'factorization_done.txt' with open(done_file, 'w') as f2: f2.write('done') @@ -557,7 +576,7 @@ cdef class MKLPardisoSolver: except FileExistsError: # flag file exists, wait for "done" file and read in factorization - done_file = self._flag_dir + 'factorization_done.txt' + done_file = self._flag_dir.decode("utf-8") + 'factorization_done.txt' while not os.path.isfile(done_file): time.sleep(1) diff --git a/tests/test.py b/tests/test.py index c71f05e..ae73545 100644 --- a/tests/test.py +++ b/tests/test.py @@ -110,6 +110,44 @@ def test_multiple_RHS(): assert rel_err < 1E3*eps return rel_err +def test_multiple_RHS_store_factorization(): + A = A_real_dict["real_symmetric_positive_definite"] + x = np.c_[xr, xr] + b = A @ x + + solver = Solver(A, "real_symmetric_positive_definite", store_factorization_dir='./') + x2 = solver.solve(b) + + eps = np.finfo(np.float64).eps + rel_err = np.linalg.norm(x-x2)/np.linalg.norm(x) + assert rel_err < 1E3*eps + return rel_err + +def test_multiple_RHS_store_factorization_clean_flag_files(): + A = A_real_dict["real_symmetric_positive_definite"] + x = np.c_[xr, xr] + b = A @ x + + solver = Solver(A, "real_symmetric_positive_definite", store_factorization_dir='./') + x2 = solver.solve(b) + + eps = np.finfo(np.float64).eps + rel_err = np.linalg.norm(x-x2)/np.linalg.norm(x) + + assert rel_err < 1E3*eps + + # run again to make sure the created flag files are checked and removed and running again works + x3 = solver.solve(b) + + eps3 = np.finfo(np.float64).eps + rel_err3 = np.linalg.norm(x-x2)/np.linalg.norm(x) + + assert rel_err3 < 1E3*eps3 + + assert rel_err == rel_err3 + + return rel_err + def test_matrix_type_errors(): A = A_real_dict["real_symmetric_positive_definite"] From 68ed2778b6d5b5ba7bef338dba151208c41fad6c Mon Sep 17 00:00:00 2001 From: JKutt Date: Tue, 3 May 2022 15:15:01 -0700 Subject: [PATCH 7/7] added pathlib for path creation as suggested by Santiago --- pydiso/mkl_solver.pyx | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/pydiso/mkl_solver.pyx b/pydiso/mkl_solver.pyx index 751c524..4328010 100644 --- a/pydiso/mkl_solver.pyx +++ b/pydiso/mkl_solver.pyx @@ -8,6 +8,7 @@ from time import time import numpy as np import scipy.sparse as sp import os +from pathlib import Path ctypedef long long MKL_INT64 ctypedef unsigned long long MKL_UINT64 @@ -319,11 +320,11 @@ cdef class MKLPardisoSolver: if store_factorization_dir is not None: # check if the flag files exist. If so delete them so factorization file get overwritten - check_file = store_factorization_dir + 'factorization_done.txt' + check_file = Path(store_factorization_dir) / 'factorization_done.txt' if os.path.exists(check_file): - second_file_to_remove = store_factorization_dir + "flagfile.txt" + second_file_to_remove = Path(store_factorization_dir) / "flagfile.txt" os.remove(check_file) os.remove(second_file_to_remove)