pardiso_wrapper/
lib.rs

1#![allow(clippy::upper_case_acronyms)]
2
3//! # PARDISO Wrapper for Rust
4//!
5//! This crate dynamically loads the PARDISO sparse solver library and provides a safe
6//! Rust interface.  It supports either MKL or Panua Pardiso backends through feature flags:
7//!
8//! - `mkl`: Intel MKL implementation (x86_64 only)
9//! - `panua`: Panua implementation
10//!
11//! Both options are supported via the common [`PardisoInterface`] trait.
12//!
13//! ### MKL Pardiso
14//!
15//! To enable dynamic linking to [`MKL Pardiso`](https://www.intel.com/content/www/us/en/docs/onemkl/developer-reference-c/),
16//! the MKL Pardiso libary (e.g. `libmkl_rt.so`) must be on the system library path
17//! (e.g. on `LD_LIBRARY_PATH` on Linux).    Alternatively, set the `MKLROOT` environment
18//! variable to the root of the MKL installation or `MKL_PARDISO_PATH` to the location
19//! of the library.  
20//!
21//! ### Panua Pardiso
22//!
23//! To enable dynamic linking to [`Panua Pardiso`](https://panua.ch/pardiso/),
24//! the Panua Pardiso library (e.g. `libpardiso.so`) must be on the system library path
25//! (e.g. on `LD_LIBRARY_PATH` on Linux).  Alternatively, set the `PARDISO_PATH` environment
26//! variable to the location of the library.
27//!
28//! Panua Pardiso is a commercial solver and requires a separate license.
29//!
30//! ## Example
31//! ```rust, no_test
32#![doc = include_str!("../examples/symmetric.rs")]
33//! ```
34
35mod enums;
36pub use enums::*;
37mod error_types;
38pub use error_types::*;
39
40cfg_if::cfg_if! {
41    if #[cfg(feature = "panua")]{
42        mod panua;
43        pub use panua::PanuaPardisoSolver;
44}}
45
46cfg_if::cfg_if! {
47    if #[cfg(feature = "mkl")]{
48        mod mkl;
49        pub use mkl::MKLPardisoSolver;
50}}
51
52#[cfg(test)]
53mod tests;
54
55#[derive(Debug)]
56#[doc(hidden)]
57pub struct PardisoData {
58    pub pt: [isize; 64],
59    pub iparm: [i32; 64],
60    pub mtype: MatrixType,
61    pub solver: SolverType,
62    pub phase: Phase,
63    pub msglvl: MessageLevel,
64    pub maxfct: i32,
65    pub mnum: i32,
66    pub perm: Vec<i32>,
67}
68
69impl Default for PardisoData {
70    fn default() -> Self {
71        Self {
72            pt: [0; 64],
73            iparm: [0; 64],
74            mtype: MatrixType::default(),
75            solver: SolverType::default(),
76            phase: Phase::default(),
77            msglvl: MessageLevel::default(),
78            maxfct: 1,
79            mnum: 1,
80            perm: vec![],
81        }
82    }
83}
84
85pub trait PardisoInterface {
86    // Getters and Setters (default implementations)
87
88    fn get_matrix_type(&self) -> MatrixType {
89        self.data().mtype
90    }
91    fn set_matrix_type(&mut self, mtype: MatrixType) {
92        self.data_mut().mtype = mtype;
93    }
94    fn get_solver(&self) -> SolverType {
95        self.data().solver
96    }
97    fn set_solver(&mut self, solver: SolverType) {
98        self.data_mut().solver = solver;
99    }
100    fn get_phase(&self) -> Phase {
101        self.data().phase
102    }
103    fn set_phase(&mut self, phase: Phase) {
104        self.data_mut().phase = phase;
105    }
106    fn get_message_level(&self) -> MessageLevel {
107        self.data().msglvl
108    }
109    fn set_message_level(&mut self, msglvl: MessageLevel) {
110        self.data_mut().msglvl = msglvl;
111    }
112    fn get_maxfct(&self) -> i32 {
113        self.data().maxfct
114    }
115    fn set_maxfct(&mut self, maxfct: i32) {
116        self.data_mut().maxfct = maxfct;
117    }
118    fn get_mnum(&self) -> i32 {
119        self.data().mnum
120    }
121    fn set_mnum(&mut self, mnum: i32) {
122        self.data_mut().mnum = mnum;
123    }
124    fn get_perm(&self) -> &[i32] {
125        self.data().perm.as_slice()
126    }
127    // : Problem here I think.   Should not own the vector passed
128    fn set_perm(&mut self, perm: &[i32]) {
129        self.data_mut().perm.resize(perm.len(), 0_i32);
130        self.data_mut().perm.copy_from_slice(perm);
131    }
132    fn get_iparm(&self, i: usize) -> i32 {
133        self.data().iparm[i]
134    }
135    fn get_iparms(&self) -> &[i32; 64] {
136        &self.data().iparm
137    }
138    fn set_iparm(&mut self, i: usize, value: i32) {
139        self.data_mut().iparm[i] = value;
140    }
141    fn get_num_positive_eigenvalues(&self) -> i32 {
142        self.data().iparm[21]
143    }
144    fn get_num_negative_eigenvalues(&self) -> i32 {
145        self.data().iparm[22]
146    }
147
148    // NB: implementors should also implement Drop and call
149    // release() in the drop method.
150    fn release(&mut self) {
151        // Set the phase to release all resources
152        self.set_phase(crate::Phase::ReleaseAll);
153
154        // call with dummies since we are releasing resource only
155        let a: Vec<f64> = vec![];
156        let ia: Vec<i32> = vec![];
157        let ja: Vec<i32> = vec![];
158        let mut b: Vec<f64> = vec![];
159        let mut x: Vec<f64> = vec![];
160
161        // Call Pardiso to release resources, ignoring any errors
162        let _ = self.pardiso(&a, &ia, &ja, &mut b, &mut x, 0, 0);
163    }
164
165    // MKL/Panua specific functions
166    #[doc(hidden)]
167    fn data(&self) -> &PardisoData;
168    #[doc(hidden)]
169    fn data_mut(&mut self) -> &mut PardisoData;
170
171    fn name(&self) -> &'static str;
172
173    fn new() -> Result<Self, PardisoError>
174    where
175        Self: Sized;
176
177    fn pardisoinit(&mut self) -> Result<(), PardisoError>;
178
179    #[allow(clippy::too_many_arguments)]
180    fn pardiso(
181        &mut self,
182        a: &[f64],
183        ia: &[i32],
184        ja: &[i32],
185        b: &mut [f64],
186        x: &mut [f64],
187        n: i32,
188        nrhs: i32,
189    ) -> Result<(), PardisoError>;
190
191    fn is_licensed() -> bool
192    where
193        Self: Sized;
194    fn is_loaded() -> bool
195    where
196        Self: Sized;
197    fn is_available() -> bool
198    where
199        Self: Sized,
200    {
201        Self::is_licensed() && Self::is_loaded()
202    }
203
204    // NB: there is no set_num_threads in this trait since
205    // Panua does not allow for configurable thread counts
206    // other than via the OMP_NUM_THREADS environment variable
207    fn get_num_threads(&self) -> Result<i32, PardisoError>;
208}
209
210// default dylib env by platform
211#[allow(dead_code)] // if no features are set
212pub(crate) fn dylib_path_env() -> &'static str {
213    cfg_if::cfg_if! {
214        if #[cfg(target_os = "macos")] {
215            "DYLD_LIBRARY_PATH"
216        } else if #[cfg(target_os = "windows")] {
217            "PATH"
218        } else {
219            "LD_LIBRARY_PATH" // linux/unknown platform
220        }
221    }
222}