rlst_suitesparse/umfpack/
complex.rs

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

rlst_suitesparse/umfpack/
complex.rs
