use crate::internal_prelude::*;
use crate::perm::{swap_cols_idx, swap_rows_idx};
use crate::utils::thread::par_split_indices;
use faer_traits::{Real, RealReg};
use linalg::matmul::matmul;
use pulp::Simd;
#[inline(always)]
fn best_value<T: ComplexField, S: Simd>(
simd: &SimdCtx<T, S>,
best_value: RealReg<T::SimdVec<S>>,
best_indices: T::SimdIndex<S>,
value: T::SimdVec<S>,
indices: T::SimdIndex<S>,
) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
let value = simd.abs1(value);
let is_better = (**simd).gt(value, best_value);
(
RealReg(simd.select(is_better, value.0, best_value.0)),
simd.iselect(is_better, indices, best_indices),
)
}
#[inline(always)]
fn best_score<T: ComplexField, S: Simd>(
simd: &SimdCtx<T, S>,
best_score: RealReg<T::SimdVec<S>>,
best_indices: T::SimdIndex<S>,
score: RealReg<T::SimdVec<S>>,
indices: T::SimdIndex<S>,
) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
let is_better = (**simd).gt(score, best_score);
(
RealReg(simd.select(is_better, score.0, best_score.0)),
simd.iselect(is_better, indices, best_indices),
)
}
#[inline(always)]
fn best_score_2d<T: ComplexField, S: Simd>(
simd: &SimdCtx<T, S>,
best_score: RealReg<T::SimdVec<S>>,
best_row: T::SimdIndex<S>,
best_col: T::SimdIndex<S>,
score: RealReg<T::SimdVec<S>>,
row: T::SimdIndex<S>,
col: T::SimdIndex<S>,
) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>, T::SimdIndex<S>) {
let is_better = (**simd).gt(score, best_score);
(
RealReg(simd.select(is_better, score.0, best_score.0)),
simd.iselect(is_better, row, best_row),
simd.iselect(is_better, col, best_col),
)
}
#[inline(always)]
fn reduce_2d<T: ComplexField, S: Simd>(
simd: &SimdCtx<T, S>,
best_values: RealReg<T::SimdVec<S>>,
best_row: T::SimdIndex<S>,
best_col: T::SimdIndex<S>,
) -> (usize, usize, Real<T>) {
let best_val = simd.reduce_max_real(best_values);
let best_val_splat = simd.splat_real(&best_val);
let is_best = (**simd).ge(best_values, best_val_splat);
let idx = simd.first_true_mask(is_best);
let best_row = bytemuck::cast_slice::<T::SimdIndex<S>, T::Index>(
core::slice::from_ref(&best_row),
)[idx];
let best_col = bytemuck::cast_slice::<T::SimdIndex<S>, T::Index>(
core::slice::from_ref(&best_col),
)[idx];
(best_row.zx(), best_col.zx(), best_val)
}
#[inline(always)]
fn best_in_col_simd<'M, T: ComplexField, S: Simd>(
simd: SimdCtx<'M, T, S>,
data: ColRef<'_, T, Dim<'M>, ContiguousFwd>,
) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
let iota = T::simd_iota(&simd.0);
let lane_count =
core::mem::size_of::<T::SimdVec<S>>() / core::mem::size_of::<T>();
let inc = simd.isplat(T::Index::truncate(lane_count));
let indices = simd.batch_indices::<4>();
let mut best_val = [RealReg(simd.zero()); 4];
let mut best_idx = [simd.izero(); 4];
let mut idx = simd.iadd(
iota,
simd.isplat(T::Index::truncate(simd.offset().wrapping_neg())),
);
simd_iter!(for (IDX, i) in [indices; 4] {
(best_val[IDX], best_idx[IDX]) = best_value(
&simd,
best_val[IDX],
best_idx[IDX],
simd.read(data, i),
idx,
);
idx = simd.iadd(idx, inc);
});
let _ = idx;
let (best_val0, best_idx0) =
best_score(&simd, best_val[0], best_idx[0], best_val[1], best_idx[1]);
let (best_val2, best_idx2) =
best_score(&simd, best_val[2], best_idx[2], best_val[3], best_idx[3]);
best_score(&simd, best_val0, best_idx0, best_val2, best_idx2)
}
#[inline(always)]
fn update_and_best_in_col_simd<'M, T: ComplexField, S: Simd>(
simd: SimdCtx<'M, T, S>,
data: ColMut<'_, T, Dim<'M>, ContiguousFwd>,
lhs: ColRef<'_, T, Dim<'M>, ContiguousFwd>,
rhs: T,
) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
let mut data = data;
let iota = T::simd_iota(&simd.0);
let lane_count =
core::mem::size_of::<T::SimdVec<S>>() / core::mem::size_of::<T>();
let inc = simd.isplat(T::Index::truncate(lane_count));
let mut best_val = [RealReg(simd.zero()); 3];
let mut best_idx = [simd.izero(); 3];
let mut idx = simd.iadd(
iota,
simd.isplat(T::Index::truncate(simd.offset().wrapping_neg())),
);
let indices = simd.batch_indices::<3>();
let rhs = simd.splat(&-rhs);
simd_iter!(for (IDX, i) in [indices; 3] {
let mut x0 = simd.read(data.rb(), i);
let l0 = simd.read(lhs, i);
x0 = simd.mul_add(l0, rhs, x0);
(best_val[IDX], best_idx[IDX]) =
best_value(&simd, best_val[IDX], best_idx[IDX], x0, idx);
simd.write(data.rb_mut(), i, x0);
idx = simd.iadd(idx, inc);
});
let _ = idx;
let (best_val0, best_idx0) =
best_score(&simd, best_val[0], best_idx[0], best_val[1], best_idx[1]);
let (best_val2, best_idx2) = (best_val[2], best_idx[2]);
best_score(&simd, best_val0, best_idx0, best_val2, best_idx2)
}
#[inline(always)]
fn best_in_mat_simd<T: ComplexField>(
data: MatRef<'_, T, usize, usize, ContiguousFwd>,
) -> (usize, usize, Real<T>) {
struct Impl<'a, 'M, 'N, T: ComplexField> {
data: MatRef<'a, T, Dim<'M>, Dim<'N>, ContiguousFwd>,
}
impl<'a, 'M, 'N, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'N, T> {
type Output = (usize, usize, Real<T>);
#[inline(always)]
fn with_simd<S: Simd>(self, simd: S) -> Self::Output {
let Self { data } = self;
let M = data.nrows();
let N = data.ncols();
let simd = SimdCtx::<'_, T, S>::new(T::simd_ctx(simd), M);
let mut best_row = simd.isplat(T::Index::truncate(0));
let mut best_col = simd.isplat(T::Index::truncate(0));
let mut best_val = simd.splat_real(&zero());
for j in N.indices() {
let col = data.col(j);
let (best_val_j, best_row_j) = best_in_col_simd(simd, col);
(best_val, best_row, best_col) = best_score_2d(
&simd,
best_val,
best_row,
best_col,
best_val_j,
best_row_j,
simd.isplat(T::Index::truncate(*j)),
);
}
reduce_2d(&simd, best_val, best_row, best_col)
}
}
with_dim!(M, data.nrows());
with_dim!(N, data.ncols());
dispatch!(
Impl {
data: data.as_shape(M, N)
},
Impl,
T
)
}
#[inline(always)]
fn update_and_best_in_mat_simd<T: ComplexField>(
data: MatMut<'_, T, usize, usize, ContiguousFwd>,
lhs: ColRef<'_, T, usize, ContiguousFwd>,
rhs: RowRef<'_, T, usize>,
align: usize,
) -> (usize, usize, Real<T>) {
struct Impl<'a, 'M, 'N, T: ComplexField> {
data: MatMut<'a, T, Dim<'M>, Dim<'N>, ContiguousFwd>,
lhs: ColRef<'a, T, Dim<'M>, ContiguousFwd>,
rhs: RowRef<'a, T, Dim<'N>>,
align: usize,
}
impl<'a, 'M, 'N, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'N, T> {
type Output = (usize, usize, Real<T>);
#[inline(always)]
fn with_simd<S: Simd>(self, simd: S) -> Self::Output {
let Self {
data,
lhs,
rhs,
align,
} = self;
let M = data.nrows();
let N = data.ncols();
let simd =
SimdCtx::<'_, T, S>::new_align(T::simd_ctx(simd), M, align);
let mut best_row = simd.isplat(T::Index::truncate(0));
let mut best_col = simd.isplat(T::Index::truncate(0));
let mut best_val = simd.splat_real(&zero());
let mut data = data;
for j in N.indices() {
let data = data.rb_mut().col_mut(j);
let rhs = rhs[j].copy();
let (best_val_j, best_row_j) =
update_and_best_in_col_simd(simd, data, lhs, rhs);
(best_val, best_row, best_col) = best_score_2d(
&simd,
best_val,
best_row,
best_col,
best_val_j,
best_row_j,
simd.isplat(T::Index::truncate(*j)),
);
}
reduce_2d(&simd, best_val, best_row, best_col)
}
}
with_dim!(M, data.nrows());
with_dim!(N, data.ncols());
dispatch!(
Impl {
data: data.as_shape_mut(M, N),
lhs: lhs.as_row_shape(M),
rhs: rhs.as_col_shape(N),
align,
},
Impl,
T
)
}
fn best_in_matrix_fallback<T: ComplexField>(
data: MatRef<'_, T>,
) -> (usize, usize, Real<T>) {
let mut max = zero();
let mut row = 0;
let mut col = 0;
let (m, n) = data.shape();
for j in 0..n {
for i in 0..m {
let abs = data[(i, j)].abs1();
if abs > max {
row = i;
col = j;
max = abs;
}
}
}
(row, col, max)
}
fn best_in_matrix<T: ComplexField>(
data: MatRef<'_, T>,
) -> (usize, usize, Real<T>) {
if const { T::SIMD_CAPABILITIES.is_simd() } {
if let Some(dst) = data.try_as_col_major() {
best_in_mat_simd(dst)
} else {
best_in_matrix_fallback(data)
}
} else {
best_in_matrix_fallback(data)
}
}
fn rank_one_update_and_best_in_matrix<T: ComplexField>(
mut dst: MatMut<'_, T>,
lhs: ColRef<'_, T>,
rhs: RowRef<'_, T>,
align: usize,
) -> (usize, usize, Real<T>) {
if const { T::SIMD_CAPABILITIES.is_simd() } {
if let (Some(dst), Some(lhs)) =
(dst.rb_mut().try_as_col_major_mut(), lhs.try_as_col_major())
{
update_and_best_in_mat_simd(dst, lhs, rhs, align)
} else {
matmul(
dst.rb_mut(),
Accum::Add,
lhs.as_mat(),
rhs.as_mat(),
-one::<T>(),
Par::Seq,
);
best_in_matrix(dst.rb())
}
} else {
matmul(
dst.rb_mut(),
Accum::Add,
lhs.as_mat(),
rhs.as_mat(),
-one::<T>(),
Par::Seq,
);
best_in_matrix(dst.rb())
}
}
fn lu_in_place_unblocked<T: ComplexField>(
A: MatMut<'_, T>,
row_trans: &mut [usize],
col_trans: &mut [usize],
par: Par,
transpose: bool,
params: Spec<FullPivLuParams, T>,
) -> usize {
let params = params.config;
let mut n_trans = 0;
let (m, n) = A.shape();
if m == 0 || n == 0 {
return 0;
}
let mut par = par;
let mut A = A;
let (mut max_row, mut max_col, mut max_score) = best_in_matrix(A.rb());
for k in 0..Ord::min(m, n) {
if max_score < min_positive() {
for (i, (row, col)) in
core::iter::zip(&mut row_trans[k..], &mut col_trans[k..])
.enumerate()
{
*row = i + k;
*col = i + k;
}
break;
}
row_trans[k] = max_row;
col_trans[k] = max_col;
if max_row != k {
swap_rows_idx(A.rb_mut(), k, max_row);
n_trans += 1;
}
if max_col != k {
swap_cols_idx(A.rb_mut(), k, max_col);
n_trans += 1;
}
let ref inv = A[(k, k)].recip();
if transpose {
for j in k + 1..n {
A[(k, j)] *= inv;
}
} else {
for i in k + 1..m {
A[(i, k)] *= inv;
}
}
if k + 1 == Ord::min(m, n) {
break;
}
if (m - k - 1) * (n - k - 1) < params.par_threshold {
par = Par::Seq;
}
let (_, A01, A10, mut A11) = A.rb_mut().split_at_mut(k + 1, k + 1);
let lhs = A10.col(k);
let rhs = A01.row(k);
match par {
Par::Seq => {
(max_row, max_col, max_score) =
rank_one_update_and_best_in_matrix(
A11.rb_mut(),
lhs,
rhs,
simd_align(k + 1),
);
},
#[cfg(feature = "rayon")]
Par::Rayon(nthreads) => {
use rayon::prelude::*;
let nthreads = nthreads.get();
let mut best = core::iter::repeat_with(|| (0, 0, zero()))
.take(nthreads)
.collect::<alloc::vec::Vec<_>>();
let full_cols = A11.ncols();
spindle::for_each(
nthreads,
best.par_iter_mut()
.zip_eq(A11.rb_mut().par_col_partition_mut(nthreads))
.zip_eq(rhs.par_partition(nthreads))
.enumerate(),
|(idx, (((max_row, max_col, max_score), A11), rhs))| {
(*max_row, *max_col, *max_score) = {
let (a, mut b, c) =
rank_one_update_and_best_in_matrix(
A11,
lhs,
rhs,
simd_align(k + 1),
);
b += par_split_indices(full_cols, idx, nthreads).0;
(a, b, c)
};
},
);
max_row = 0;
max_col = 0;
max_score = zero();
for (row, col, val) in best {
if val > max_score {
max_row = row;
max_col = col;
max_score = val;
}
}
},
}
max_row += k + 1;
max_col += k + 1;
}
n_trans
}
#[derive(Copy, Clone, Debug)]
pub struct FullPivLuParams {
pub par_threshold: usize,
#[doc(hidden)]
pub non_exhaustive: NonExhaustive,
}
impl<T: ComplexField> Auto<T> for FullPivLuParams {
#[inline]
fn auto() -> Self {
Self {
par_threshold: 256 * 512,
non_exhaustive: NonExhaustive(()),
}
}
}
#[inline]
pub fn lu_in_place_scratch<I: Index, T: ComplexField>(
nrows: usize,
ncols: usize,
par: Par,
params: Spec<FullPivLuParams, T>,
) -> StackReq {
_ = par;
_ = params;
let size = Ord::min(nrows, ncols);
StackReq::new::<usize>(size).array(2)
}
#[derive(Copy, Clone, Debug)]
pub struct FullPivLuInfo {
pub transposition_count: usize,
}
pub fn lu_in_place<'out, I: Index, T: ComplexField>(
mat: MatMut<'_, T>,
row_perm: &'out mut [I],
row_perm_inv: &'out mut [I],
col_perm: &'out mut [I],
col_perm_inv: &'out mut [I],
par: Par,
stack: &mut MemStack,
params: Spec<FullPivLuParams, T>,
) -> (FullPivLuInfo, PermRef<'out, I>, PermRef<'out, I>) {
#[cfg(feature = "perf-warn")]
if (mat.col_stride().unsigned_abs() == 1
|| mat.row_stride().unsigned_abs() != 1)
&& crate::__perf_warn!(LU_WARN)
{
log::warn!(
target : "faer_perf",
"LU with full pivoting prefers column-major or row-major matrix. Found matrix with generic strides."
);
}
let (M, N) = mat.shape();
let size = Ord::min(M, N);
let (mut row_transpositions, stack) = stack.make_with(size, |_| 0);
let row_transpositions = row_transpositions.as_mut();
let (mut col_transpositions, _) = stack.make_with(size, |_| 0);
let col_transpositions = col_transpositions.as_mut();
let n_transpositions = if mat.row_stride().abs() < mat.col_stride().abs() {
lu_in_place_unblocked(
mat,
row_transpositions,
col_transpositions,
par,
false,
params,
)
} else {
lu_in_place_unblocked(
mat.transpose_mut(),
col_transpositions,
row_transpositions,
par,
true,
params,
)
};
for i in 0..M {
row_perm[i] = I::truncate(i);
}
for (i, t) in row_transpositions.iter().copied().enumerate() {
row_perm.as_mut().swap(i, t);
}
for i in 0..M {
row_perm_inv[row_perm[i].zx()] = I::truncate(i);
}
for j in 0..N {
col_perm[j] = I::truncate(j);
}
for (i, t) in col_transpositions.iter().copied().enumerate() {
col_perm.as_mut().swap(i, t);
}
for j in 0..N {
col_perm_inv[col_perm[j].zx()] = I::truncate(j);
}
unsafe {
(
FullPivLuInfo {
transposition_count: n_transpositions,
},
PermRef::new_unchecked(row_perm, row_perm_inv, M),
PermRef::new_unchecked(col_perm, col_perm_inv, N),
)
}
}
#[cfg(test)]
mod tests {
use super::*;
use crate::stats::prelude::*;
use crate::utils::approx::*;
use crate::{Mat, assert, c64};
use dyn_stack::MemBuffer;
#[test]
fn test_flu() {
let rng = &mut StdRng::seed_from_u64(0);
for par in [Par::Seq, Par::rayon(8)] {
for m in [8, 16, 24, 32, 128, 255, 256, 257] {
let n = 8;
let approx_eq = CwiseMat(ApproxEq {
abs_tol: 1e-10,
rel_tol: 1e-10,
});
let A = CwiseMatDistribution {
nrows: m,
ncols: n,
dist: ComplexDistribution::new(
StandardNormal,
StandardNormal,
),
}
.rand::<Mat<c64>>(rng);
let A = A.as_ref();
let mut LU = A.cloned();
let row_perm = &mut *vec![0usize; m];
let row_perm_inv = &mut *vec![0usize; m];
let col_perm = &mut *vec![0usize; n];
let col_perm_inv = &mut *vec![0usize; n];
let (_, p, q) = lu_in_place(
LU.as_mut(),
row_perm,
row_perm_inv,
col_perm,
col_perm_inv,
par,
MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
usize,
c64,
>(
n, n, par, default()
))),
default(),
);
let mut L = LU.as_ref().cloned();
let mut U = LU.as_ref().cloned();
for j in 0..n {
for i in 0..j {
L[(i, j)] = c64::ZERO;
}
L[(j, j)] = c64::ONE;
}
for j in 0..n {
for i in j + 1..m {
U[(i, j)] = c64::ZERO;
}
}
let L = L.as_ref();
let U = U.as_ref();
let U = U.subrows(0, n);
assert!(p.inverse() * L * U * q ~ A);
}
for n in [16, 24, 32, 128, 255, 256, 257] {
let approx_eq = CwiseMat(ApproxEq {
abs_tol: 1e-10,
rel_tol: 1e-10,
});
let A = CwiseMatDistribution {
nrows: n,
ncols: n,
dist: StandardNormal,
}
.rand::<Mat<f64>>(rng);
let A = A.as_ref();
let mut LU = A.cloned();
let row_perm = &mut *vec![0usize; n];
let row_perm_inv = &mut *vec![0usize; n];
let col_perm = &mut *vec![0usize; n];
let col_perm_inv = &mut *vec![0usize; n];
let (_, p, q) = lu_in_place(
LU.as_mut(),
row_perm,
row_perm_inv,
col_perm,
col_perm_inv,
par,
MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
usize,
f64,
>(
n, n, par, default()
))),
default(),
);
let mut L = LU.as_ref().cloned();
let mut U = LU.as_ref().cloned();
for j in 0..n {
for i in 0..j {
L[(i, j)] = 0.0;
}
L[(j, j)] = 1.0;
}
for j in 0..n {
for i in j + 1..n {
U[(i, j)] = 0.0;
}
}
let L = L.as_ref();
let U = U.as_ref();
assert!(p.inverse() * L * U * q ~ A);
}
}
}
#[test]
fn test_gh238() {
#[rustfmt::skip]
let A = [
[
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-2.0,
0.0,
1.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
1.0,
],
[
0.0,
0.0,
0.0,
1.0,
0.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
-1.0,
0.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
1.0,
0.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-1.0,
0.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
1.0,
0.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
1.0,
0.0,
0.0,
-1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
1.0,
0.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
1.0,
0.0,
0.0,
-1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
1.0,
0.0,
0.0,
1.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
-1.0,
0.0,
0.0,
1.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
1.0,
0.0,
0.0,
1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
[
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
1.0,
0.0,
0.0,
-1.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
0.0,
],
];
let A = MatRef::from_row_major_array(&A).to_owned();
let b = [
-0.0, -0.0, -2.0, -2.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0,
-0.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
0.0, 0.0, 0.0, 0.0,
];
let b = ColRef::from_slice(&b);
let lu = A.full_piv_lu();
assert!(lu.solve(b).is_all_finite());
}
}