1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
//! Implementation of ?getri - Matrix inversion using LU factorization.
use lapack::{cgetri, dgetri, sgetri, zgetri};
use crate::base_types::LapackResult;
use crate::base_types::{LapackError, c32, c64};
use crate::dense::linalg::lapack::interface::lapack_return;
use num::{Zero, complex::ComplexFloat};
/// ?getri - Matrix inversion using LU factorization.
///
/// **Arguments:**
///
/// - `n`: Number of rows and columns in the square matrix `a`.
/// - `a`: The factors `L` and `U` from the LU factorization of the matrix.
/// - `ipiv`: Pivot indices from the LU factorization.
/// - `lda`: Leading dimension of the matrix `a`.
///
/// **Returns:**
/// - The Lapack info code if the inversion fails.
pub trait Getri: Sized {
/// Compute the inverse of a matrix `a` using the LU factors and pivot indices.
fn getri(n: usize, a: &mut [Self], lda: usize, ipiv: &[i32]) -> LapackResult<()>;
}
macro_rules! implement_getri {
($scalar:ty, $getri:expr) => {
impl Getri for $scalar {
fn getri(n: usize, a: &mut [Self], lda: usize, ipiv: &[i32]) -> LapackResult<()> {
let mut info = 0;
assert!(
lda >= std::cmp::max(1, n),
"Require `lda` {} >= `max(1, N)` {} .",
lda,
std::cmp::max(1, n)
);
assert_eq!(
a.len(),
lda * n,
"Require `a.len()` {} == `lda * n` {}.",
a.len(),
lda * n,
);
assert_eq!(
ipiv.len(),
n,
"Require `ipiv.len()` {} == `n` {}.",
ipiv.len(),
n,
);
let mut work = vec![<$scalar>::zero(); 1];
unsafe {
$getri(n as i32, a, lda as i32, ipiv, &mut work, -1, &mut info);
}
let lwork = work[0].re() as i32;
if info != 0 {
return Err(LapackError::LapackInfoCode(info));
}
let mut work = vec![<$scalar>::zero(); lwork as usize];
unsafe {
$getri(n as i32, a, lda as i32, ipiv, &mut work, lwork, &mut info);
}
lapack_return(info, ())
}
}
};
}
implement_getri!(f32, sgetri);
implement_getri!(f64, dgetri);
implement_getri!(c32, cgetri);
implement_getri!(c64, zgetri);