mdarray_linalg_lapack/solve/
context.rs

1//! Linear System Solver using LU decomposition (GESV):
2//!     AX = B
3//! where:
4//!     - A is n × n (square coefficient matrix, overwritten with LU factorization)
5//!     - X is n × nrhs (solution matrix)
6//!     - B is n × nrhs (right-hand side matrix, overwritten with solution)
7//!
8//! The function `gesv` (LAPACK) solves a system of linear equations AX = B using LU decomposition with partial pivoting.
9//! It computes the LU factorization of A and then uses it to solve the linear system.
10//! The matrix A is overwritten by its LU factorization, and B is overwritten by the solution X.
11
12use mdarray::{Array, Dense, Dim, Layout, Shape, Slice};
13use mdarray_linalg::solve::{Solve, SolveError};
14use num_complex::ComplexFloat;
15
16use super::{scalar::LapackScalar, simple::gesv};
17use crate::Lapack;
18
19impl<T, D: Dim> Solve<T, D> for Lapack
20where
21    T: ComplexFloat + Default + LapackScalar,
22    T::Real: Into<T>,
23{
24    fn solve_write<R: Dim, La: Layout, Lb: Layout>(
25        &self,
26        a: &mut Slice<T, (D, D), La>,
27        b: &mut Slice<T, (D, R), Lb>,
28    ) -> Result<(), SolveError> {
29        gesv::<_, Lb, T, D, R>(a, b)
30    }
31
32    fn solve<R: Dim, La: Layout, Lb: Layout>(
33        &self,
34        a: &mut Slice<T, (D, D), La>,
35        b: &Slice<T, (D, R), Lb>,
36    ) -> Result<Array<T, (D, R)>, SolveError> {
37        let ash = *a.shape();
38        let bsh = *b.shape();
39
40        let n = ash.dim(0);
41        let nrhs = bsh.dim(1);
42
43        let mut b_copy = Array::from_elem(<(D, R) as Shape>::from_dims(&[n, nrhs]), T::default());
44
45        for i in 0..n {
46            for j in 0..nrhs {
47                b_copy[[i, j]] = b[[i, j]];
48            }
49        }
50
51        gesv::<_, Dense, T, D, R>(a, &mut b_copy)?;
52        Ok(b_copy)
53    }
54}
mdarray_linalg_lapack/solve/context.rs

mdarray_linalg_lapack/solve/
context.rs
