rlst_suitesparse/umfpack/
real.rs

1use std::ffi::c_void;
2use std::ptr;
3use std::sync::Arc;
4
5use rlst::sparse::csr_mat::CsrMatrix;
6use rlst::{Array, DynArray, MutableRawAccessArrayImpl, RawAccessArrayImpl, rlst_dynamic_array};
7
8use super::csc::CscMatrix;
9use super::sys;
10use super::{
11    UMFPACK_A, UmfpackError, UmfpackFactorizationStatus, UmfpackFactorize, UmfpackResult,
12    array_data, array_data_mut, factorization_status,
13};
14
15struct UmfpackRealSymbolic {
16    raw: *mut c_void,
17    matrix: Arc<CscMatrix<f64>>,
18}
19
20impl UmfpackFactorize for CsrMatrix<f64> {
21    type Numeric = UmfpackNumeric;
22
23    fn factorize(&self) -> UmfpackResult<Self::Numeric> {
24        UmfpackRealSymbolic::analyze(self)?.factorize()
25    }
26}
27
28impl UmfpackRealSymbolic {
29    fn analyze(matrix: &CsrMatrix<f64>) -> UmfpackResult<Self> {
30        let matrix = Arc::new(CscMatrix::from_csr(matrix)?);
31        let mut raw = ptr::null_mut();
32        let status = unsafe {
33            sys::umfpack_dl_symbolic(
34                matrix.rows_i64()?,
35                matrix.columns_i64()?,
36                matrix.column_pointers.as_ptr(),
37                matrix.row_indices.as_ptr(),
38                matrix.values.as_ptr(),
39                &mut raw,
40                ptr::null(),
41                ptr::null_mut(),
42            )
43        };
44
45        if status != 0 {
46            if !raw.is_null() {
47                unsafe { sys::umfpack_dl_free_symbolic(&mut raw) };
48            }
49            return Err(UmfpackError::NativeFailure {
50                operation: "symbolic analysis",
51                status,
52            });
53        }
54        if raw.is_null() {
55            return Err(UmfpackError::NullHandle("symbolic analysis"));
56        }
57
58        Ok(Self { raw, matrix })
59    }
60
61    fn factorize(&self) -> UmfpackResult<UmfpackNumeric> {
62        let mut raw = ptr::null_mut();
63        let status = unsafe {
64            sys::umfpack_dl_numeric(
65                self.matrix.column_pointers.as_ptr(),
66                self.matrix.row_indices.as_ptr(),
67                self.matrix.values.as_ptr(),
68                self.raw,
69                &mut raw,
70                ptr::null(),
71                ptr::null_mut(),
72            )
73        };
74
75        if status < 0 {
76            if !raw.is_null() {
77                unsafe { sys::umfpack_dl_free_numeric(&mut raw) };
78            }
79            return Err(UmfpackError::NativeFailure {
80                operation: "numeric factorization",
81                status,
82            });
83        }
84        if raw.is_null() {
85            return Err(UmfpackError::NullHandle("numeric factorization"));
86        }
87
88        Ok(UmfpackNumeric {
89            raw,
90            matrix: Arc::clone(&self.matrix),
91            status: factorization_status(status),
92        })
93    }
94}
95
96impl Drop for UmfpackRealSymbolic {
97    fn drop(&mut self) {
98        if !self.raw.is_null() {
99            unsafe { sys::umfpack_dl_free_symbolic(&mut self.raw) };
100        }
101    }
102}
103
104/// Numeric UMFPACK factorization of an RLST CSR matrix.
105pub struct UmfpackNumeric {
106    raw: *mut c_void,
107    matrix: Arc<CscMatrix<f64>>,
108    status: UmfpackFactorizationStatus,
109}
110
111impl UmfpackNumeric {
112    /// Return the non-error status reported while this factorization was
113    /// created.
114    pub fn status(&self) -> UmfpackFactorizationStatus {
115        self.status
116    }
117
118    /// Solve `Ax = rhs` and return the solution in a new contiguous dynamic array.
119    ///
120    /// `rhs` must expose contiguous storage and have length equal to the matrix order.
121    pub fn solve<ArrayImpl>(&self, rhs: &Array<ArrayImpl, 1>) -> UmfpackResult<DynArray<f64, 1>>
122    where
123        ArrayImpl: RawAccessArrayImpl<f64, 1>,
124    {
125        let mut solution = rlst_dynamic_array!(f64, [self.matrix.rows]);
126        self.solve_into(rhs, &mut solution)?;
127        Ok(solution)
128    }
129
130    /// Solve `Ax = rhs`, writing the result into `solution`.
131    ///
132    /// Both arrays must expose contiguous storage and have length equal to the
133    /// matrix order. Native solve warnings and failures are returned as
134    /// [`UmfpackError::NativeFailure`].
135    pub fn solve_into<RhsArrayImpl, SolutionArrayImpl>(
136        &self,
137        rhs: &Array<RhsArrayImpl, 1>,
138        solution: &mut Array<SolutionArrayImpl, 1>,
139    ) -> UmfpackResult<()>
140    where
141        RhsArrayImpl: RawAccessArrayImpl<f64, 1>,
142        SolutionArrayImpl: MutableRawAccessArrayImpl<f64, 1>,
143    {
144        let rhs = array_data(rhs, "right-hand side")?;
145        let solution = array_data_mut(solution, "solution")?;
146        self.validate_vector_len(rhs.len())?;
147        self.validate_vector_len(solution.len())?;
148
149        let status = unsafe {
150            sys::umfpack_dl_solve(
151                UMFPACK_A,
152                self.matrix.column_pointers.as_ptr(),
153                self.matrix.row_indices.as_ptr(),
154                self.matrix.values.as_ptr(),
155                solution.as_mut_ptr(),
156                rhs.as_ptr(),
157                self.raw,
158                ptr::null(),
159                ptr::null_mut(),
160            )
161        };
162
163        if status == 0 {
164            Ok(())
165        } else {
166            Err(UmfpackError::NativeFailure {
167                operation: "solve",
168                status,
169            })
170        }
171    }
172
173    fn validate_vector_len(&self, actual: usize) -> UmfpackResult<()> {
174        if actual == self.matrix.rows {
175            Ok(())
176        } else {
177            Err(UmfpackError::DimensionMismatch {
178                expected: self.matrix.rows,
179                actual,
180            })
181        }
182    }
183}
184
185impl Drop for UmfpackNumeric {
186    fn drop(&mut self) {
187        if !self.raw.is_null() {
188            unsafe { sys::umfpack_dl_free_numeric(&mut self.raw) };
189        }
190    }
191}
rlst_suitesparse/umfpack/real.rs

rlst_suitesparse/umfpack/
real.rs
