1use crate::internal_prelude::*;
2use crate::utils::cost::{ParThreshold, rate};
3use crate::perm::{swap_cols_idx, swap_rows_idx};
4use crate::utils::thread::par_split_indices;
5use faer_traits::{Real, RealReg};
6use linalg::matmul::matmul;
7use pulp::Simd;
8
9#[inline(always)]
10fn best_value<T: ComplexField, S: Simd>(
11 simd: &SimdCtx<T, S>,
12 best_value: RealReg<T::SimdVec<S>>,
13 best_indices: T::SimdIndex<S>,
14 value: T::SimdVec<S>,
15 indices: T::SimdIndex<S>,
16) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
17 let value = simd.abs1(value);
18 let is_better = (**simd).gt(value, best_value);
19 (
20 RealReg(simd.select(is_better, value.0, best_value.0)),
21 simd.iselect(is_better, indices, best_indices),
22 )
23}
24#[inline(always)]
25fn best_score<T: ComplexField, S: Simd>(
26 simd: &SimdCtx<T, S>,
27 best_score: RealReg<T::SimdVec<S>>,
28 best_indices: T::SimdIndex<S>,
29 score: RealReg<T::SimdVec<S>>,
30 indices: T::SimdIndex<S>,
31) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
32 let is_better = (**simd).gt(score, best_score);
33 (
34 RealReg(simd.select(is_better, score.0, best_score.0)),
35 simd.iselect(is_better, indices, best_indices),
36 )
37}
38#[inline(always)]
39fn best_score_2d<T: ComplexField, S: Simd>(
40 simd: &SimdCtx<T, S>,
41 best_score: RealReg<T::SimdVec<S>>,
42 best_row: T::SimdIndex<S>,
43 best_col: T::SimdIndex<S>,
44 score: RealReg<T::SimdVec<S>>,
45 row: T::SimdIndex<S>,
46 col: T::SimdIndex<S>,
47) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>, T::SimdIndex<S>) {
48 let is_better = (**simd).gt(score, best_score);
49 (
50 RealReg(simd.select(is_better, score.0, best_score.0)),
51 simd.iselect(is_better, row, best_row),
52 simd.iselect(is_better, col, best_col),
53 )
54}
55#[inline(always)]
56fn reduce_2d<T: ComplexField, S: Simd>(
57 simd: &SimdCtx<T, S>,
58 best_values: RealReg<T::SimdVec<S>>,
59 best_row: T::SimdIndex<S>,
60 best_col: T::SimdIndex<S>,
61) -> (usize, usize, Real<T>) {
62 let best_val = simd.reduce_max_real(best_values);
63 let best_val_splat = simd.splat_real(&best_val);
64 let is_best = (**simd).ge(best_values, best_val_splat);
65 let idx = simd.first_true_mask(is_best);
66 let best_row = bytemuck::cast_slice::<T::SimdIndex<S>, T::Index>(
67 core::slice::from_ref(&best_row),
68 )[idx];
69 let best_col = bytemuck::cast_slice::<T::SimdIndex<S>, T::Index>(
70 core::slice::from_ref(&best_col),
71 )[idx];
72 (best_row.zx(), best_col.zx(), best_val)
73}
74#[inline(always)]
75fn best_in_col_simd<'M, T: ComplexField, S: Simd>(
76 simd: SimdCtx<'M, T, S>,
77 data: ColRef<'_, T, Dim<'M>, ContiguousFwd>,
78) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
79 let iota = T::simd_iota(&simd.0);
80 let lane_count =
81 core::mem::size_of::<T::SimdVec<S>>() / core::mem::size_of::<T>();
82 let inc = simd.isplat(T::Index::truncate(lane_count));
83
84 let indices = simd.batch_indices::<4>();
85 let mut best_val = [RealReg(simd.zero()); 4];
86 let mut best_idx = [simd.izero(); 4];
87
88 let mut idx = simd.iadd(
89 iota,
90 simd.isplat(T::Index::truncate(simd.offset().wrapping_neg())),
91 );
92 simd_iter!(for (IDX, i) in [indices; 4] {
93 (best_val[IDX], best_idx[IDX]) = best_value(
94 &simd,
95 best_val[IDX],
96 best_idx[IDX],
97 simd.read(data, i),
98 idx,
99 );
100 idx = simd.iadd(idx, inc);
101 });
102 let _ = idx;
103
104 let (best_val0, best_idx0) =
105 best_score(&simd, best_val[0], best_idx[0], best_val[1], best_idx[1]);
106 let (best_val2, best_idx2) =
107 best_score(&simd, best_val[2], best_idx[2], best_val[3], best_idx[3]);
108 best_score(&simd, best_val0, best_idx0, best_val2, best_idx2)
109}
110#[inline(always)]
111fn update_and_best_in_col_simd<'M, T: ComplexField, S: Simd>(
112 simd: SimdCtx<'M, T, S>,
113 data: ColMut<'_, T, Dim<'M>, ContiguousFwd>,
114 lhs: ColRef<'_, T, Dim<'M>, ContiguousFwd>,
115 rhs: T,
116) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
117 let mut data = data;
118 let iota = T::simd_iota(&simd.0);
119 let lane_count =
120 core::mem::size_of::<T::SimdVec<S>>() / core::mem::size_of::<T>();
121 let inc = simd.isplat(T::Index::truncate(lane_count));
122 let mut best_val = [RealReg(simd.zero()); 3];
123 let mut best_idx = [simd.izero(); 3];
124
125 let mut idx = simd.iadd(
126 iota,
127 simd.isplat(T::Index::truncate(simd.offset().wrapping_neg())),
128 );
129
130 let indices = simd.batch_indices::<3>();
131 let rhs = simd.splat(&-rhs);
132 simd_iter!(for (IDX, i) in [indices; 3] {
133 let mut x0 = simd.read(data.rb(), i);
134 let l0 = simd.read(lhs, i);
135 x0 = simd.mul_add(l0, rhs, x0);
136 (best_val[IDX], best_idx[IDX]) =
137 best_value(&simd, best_val[IDX], best_idx[IDX], x0, idx);
138 simd.write(data.rb_mut(), i, x0);
139 idx = simd.iadd(idx, inc);
140 });
141 let _ = idx;
142
143 let (best_val0, best_idx0) =
144 best_score(&simd, best_val[0], best_idx[0], best_val[1], best_idx[1]);
145 let (best_val2, best_idx2) = (best_val[2], best_idx[2]);
146 best_score(&simd, best_val0, best_idx0, best_val2, best_idx2)
147}
148#[inline(always)]
149fn best_in_mat_simd<T: ComplexField>(
150 data: MatRef<'_, T, usize, usize, ContiguousFwd>,
151) -> (usize, usize, Real<T>) {
152 struct Impl<'a, 'M, 'N, T: ComplexField> {
153 data: MatRef<'a, T, Dim<'M>, Dim<'N>, ContiguousFwd>,
154 }
155 impl<'a, 'M, 'N, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'N, T> {
156 type Output = (usize, usize, Real<T>);
157
158 #[inline(always)]
159 fn with_simd<S: Simd>(self, simd: S) -> Self::Output {
160 let Self { data } = self;
161 let M = data.nrows();
162 let N = data.ncols();
163 let simd = SimdCtx::<'_, T, S>::new(T::simd_ctx(simd), M);
164 let mut best_row = simd.isplat(T::Index::truncate(0));
165 let mut best_col = simd.isplat(T::Index::truncate(0));
166 let mut best_val = simd.splat_real(&zero());
167 for j in N.indices() {
168 let col = data.col(j);
169 let (best_val_j, best_row_j) = best_in_col_simd(simd, col);
170 (best_val, best_row, best_col) = best_score_2d(
171 &simd,
172 best_val,
173 best_row,
174 best_col,
175 best_val_j,
176 best_row_j,
177 simd.isplat(T::Index::truncate(*j)),
178 );
179 }
180 reduce_2d(&simd, best_val, best_row, best_col)
181 }
182 }
183 with_dim!(M, data.nrows());
184 with_dim!(N, data.ncols());
185 dispatch!(
186 Impl {
187 data: data.as_shape(M, N)
188 },
189 Impl,
190 T
191 )
192}
193#[inline(always)]
194fn update_and_best_in_mat_simd<T: ComplexField>(
195 data: MatMut<'_, T, usize, usize, ContiguousFwd>,
196 lhs: ColRef<'_, T, usize, ContiguousFwd>,
197 rhs: RowRef<'_, T, usize>,
198 align: usize,
199) -> (usize, usize, Real<T>) {
200 struct Impl<'a, 'M, 'N, T: ComplexField> {
201 data: MatMut<'a, T, Dim<'M>, Dim<'N>, ContiguousFwd>,
202 lhs: ColRef<'a, T, Dim<'M>, ContiguousFwd>,
203 rhs: RowRef<'a, T, Dim<'N>>,
204 align: usize,
205 }
206 impl<'a, 'M, 'N, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'N, T> {
207 type Output = (usize, usize, Real<T>);
208
209 #[inline(always)]
210 fn with_simd<S: Simd>(self, simd: S) -> Self::Output {
211 let Self {
212 data,
213 lhs,
214 rhs,
215 align,
216 } = self;
217 let M = data.nrows();
218 let N = data.ncols();
219 let simd =
220 SimdCtx::<'_, T, S>::new_align(T::simd_ctx(simd), M, align);
221 let mut best_row = simd.isplat(T::Index::truncate(0));
222 let mut best_col = simd.isplat(T::Index::truncate(0));
223 let mut best_val = simd.splat_real(&zero());
224 let mut data = data;
225 for j in N.indices() {
226 let data = data.rb_mut().col_mut(j);
227 let rhs = rhs[j].copy();
228 let (best_val_j, best_row_j) =
229 update_and_best_in_col_simd(simd, data, lhs, rhs);
230 (best_val, best_row, best_col) = best_score_2d(
231 &simd,
232 best_val,
233 best_row,
234 best_col,
235 best_val_j,
236 best_row_j,
237 simd.isplat(T::Index::truncate(*j)),
238 );
239 }
240 reduce_2d(&simd, best_val, best_row, best_col)
241 }
242 }
243 with_dim!(M, data.nrows());
244 with_dim!(N, data.ncols());
245 dispatch!(
246 Impl {
247 data: data.as_shape_mut(M, N),
248 lhs: lhs.as_row_shape(M),
249 rhs: rhs.as_col_shape(N),
250 align,
251 },
252 Impl,
253 T
254 )
255}
256fn best_in_matrix_fallback<T: ComplexField>(
257 data: MatRef<'_, T>,
258) -> (usize, usize, Real<T>) {
259 let mut max = zero();
260 let mut row = 0;
261 let mut col = 0;
262 let (m, n) = data.shape();
263 for j in 0..n {
264 for i in 0..m {
265 let abs = data[(i, j)].abs1();
266 if abs > max {
267 row = i;
268 col = j;
269 max = abs;
270 }
271 }
272 }
273 (row, col, max)
274}
275fn best_in_matrix<T: ComplexField>(
276 data: MatRef<'_, T>,
277) -> (usize, usize, Real<T>) {
278 if const { T::SIMD_CAPABILITIES.is_simd() } {
279 if let Some(dst) = data.try_as_col_major() {
280 best_in_mat_simd(dst)
281 } else {
282 best_in_matrix_fallback(data)
283 }
284 } else {
285 best_in_matrix_fallback(data)
286 }
287}
288fn rank_one_update_and_best_in_matrix<T: ComplexField>(
289 mut dst: MatMut<'_, T>,
290 lhs: ColRef<'_, T>,
291 rhs: RowRef<'_, T>,
292 align: usize,
293) -> (usize, usize, Real<T>) {
294 if const { T::SIMD_CAPABILITIES.is_simd() } {
295 if let (Some(dst), Some(lhs)) =
296 (dst.rb_mut().try_as_col_major_mut(), lhs.try_as_col_major())
297 {
298 update_and_best_in_mat_simd(dst, lhs, rhs, align)
299 } else {
300 matmul(
301 dst.rb_mut(),
302 Accum::Add,
303 lhs.as_mat(),
304 rhs.as_mat(),
305 -one::<T>(),
306 Par::Seq,
307 );
308 best_in_matrix(dst.rb())
309 }
310 } else {
311 matmul(
312 dst.rb_mut(),
313 Accum::Add,
314 lhs.as_mat(),
315 rhs.as_mat(),
316 -one::<T>(),
317 Par::Seq,
318 );
319 best_in_matrix(dst.rb())
320 }
321}
322fn lu_in_place_unblocked<T: ComplexField>(
323 A: MatMut<'_, T>,
324 row_trans: &mut [usize],
325 col_trans: &mut [usize],
326 par: Par,
327 transpose: bool,
328 params: Spec<FullPivLuParams, T>,
329) -> usize {
330 let params = params.config;
331 let mut n_trans = 0;
332 let (m, n) = A.shape();
333 if m == 0 || n == 0 {
334 return 0;
335 }
336 let mut par = par;
337 let mut A = A;
338 let (mut max_row, mut max_col, mut max_score) = best_in_matrix(A.rb());
339 for k in 0..Ord::min(m, n) {
340 if max_score < min_positive() {
341 for (i, (row, col)) in
342 core::iter::zip(&mut row_trans[k..], &mut col_trans[k..])
343 .enumerate()
344 {
345 *row = i + k;
346 *col = i + k;
347 }
348 break;
349 }
350 row_trans[k] = max_row;
351 col_trans[k] = max_col;
352 if max_row != k {
353 swap_rows_idx(A.rb_mut(), k, max_row);
354 n_trans += 1;
355 }
356 if max_col != k {
357 swap_cols_idx(A.rb_mut(), k, max_col);
358 n_trans += 1;
359 }
360 let ref inv = A[(k, k)].recip();
361 if transpose {
362 for j in k + 1..n {
363 A[(k, j)] *= inv;
364 }
365 } else {
366 for i in k + 1..m {
367 A[(i, k)] *= inv;
368 }
369 }
370 if k + 1 == Ord::min(m, n) {
371 break;
372 }
373 if !params.par_threshold.region_pays::<T>(
376 par.shape(),
377 (m - k - 1) * (n - k - 1),
378 rate::lu_full_pivot(),
379 1,
380 ) {
381 par = Par::Seq;
382 }
383 let (_, A01, A10, mut A11) = A.rb_mut().split_at_mut(k + 1, k + 1);
384 let lhs = A10.col(k);
385 let rhs = A01.row(k);
386 match par {
387 Par::Seq => {
388 (max_row, max_col, max_score) =
389 rank_one_update_and_best_in_matrix(
390 A11.rb_mut(),
391 lhs,
392 rhs,
393 simd_align(k + 1),
394 );
395 },
396 #[cfg(feature = "rayon")]
397 _ => {
398 let nthreads = core::num::NonZeroUsize::new(par.degree()).unwrap();
399 use rayon::prelude::*;
400 let nthreads = nthreads.get();
401 let mut best = core::iter::repeat_with(|| (0, 0, zero()))
402 .take(nthreads)
403 .collect::<alloc::vec::Vec<_>>();
404 let full_cols = A11.ncols();
405 crate::utils::thread::for_each(par.reborrow(),
406 nthreads,
407 best.par_iter_mut()
408 .zip_eq(A11.rb_mut().par_col_partition_mut(nthreads))
409 .zip_eq(rhs.par_partition(nthreads))
410 .enumerate(),
411 |(idx, (((max_row, max_col, max_score), A11), rhs))| {
412 (*max_row, *max_col, *max_score) = {
413 let (a, mut b, c) =
414 rank_one_update_and_best_in_matrix(
415 A11,
416 lhs,
417 rhs,
418 simd_align(k + 1),
419 );
420 b += par_split_indices(full_cols, idx, nthreads).0;
421 (a, b, c)
422 };
423 },
424 );
425 max_row = 0;
426 max_col = 0;
427 max_score = zero();
428 for (row, col, val) in best {
429 if val > max_score {
430 max_row = row;
431 max_col = col;
432 max_score = val;
433 }
434 }
435 },
436 }
437 max_row += k + 1;
438 max_col += k + 1;
439 }
440 n_trans
441}
442#[derive(Copy, Clone, Debug)]
444pub struct FullPivLuParams {
445 pub par_threshold: ParThreshold,
449 #[doc(hidden)]
450 pub non_exhaustive: NonExhaustive,
451}
452impl<T: ComplexField> Auto<T> for FullPivLuParams {
453 #[inline]
454 fn auto() -> Self {
455 Self {
456 par_threshold: ParThreshold::Derived,
457 non_exhaustive: NonExhaustive(()),
458 }
459 }
460}
461#[inline]
462pub fn lu_in_place_scratch<I: Index, T: ComplexField>(
463 nrows: usize,
464 ncols: usize,
465 par: ParShape,
466 params: Spec<FullPivLuParams, T>,
467) -> StackReq {
468 _ = par;
469 _ = params;
470 let size = Ord::min(nrows, ncols);
471 StackReq::new::<usize>(size).array(2)
472}
473#[derive(Copy, Clone, Debug)]
474pub struct FullPivLuInfo {
475 pub transposition_count: usize,
476}
477pub fn lu_in_place<'out, I: Index, T: ComplexField>(
478 mat: MatMut<'_, T>,
479 row_perm: &'out mut [I],
480 row_perm_inv: &'out mut [I],
481 col_perm: &'out mut [I],
482 col_perm_inv: &'out mut [I],
483 mut par: Par,
484 stack: &mut MemStack,
485 params: Spec<FullPivLuParams, T>,
486) -> (FullPivLuInfo, PermRef<'out, I>, PermRef<'out, I>) {
487 #[cfg(feature = "perf-warn")]
488 if (mat.col_stride().unsigned_abs() == 1
489 || mat.row_stride().unsigned_abs() != 1)
490 && crate::__perf_warn!(LU_WARN)
491 {
492 log::warn!(
493 target : "faer_perf",
494 "LU with full pivoting prefers column-major or row-major matrix. Found matrix with generic strides."
495 );
496 }
497 let (M, N) = mat.shape();
499 if mat.row_stride().unsigned_abs() != 1
500 && M >= N
501 && N > 1
502 && linalg::row_major_copy_pays(M, N)
503 {
504 let mut work = alloc::vec![zero::<T>(); M * N];
505 let mut work = MatMut::from_column_major_slice_mut(&mut work, M, N);
506 linalg::copy_across_layouts(work.rb_mut(), mat.rb(), par.reborrow());
507 let ret = lu_in_place(
508 work.rb_mut(),
509 row_perm,
510 row_perm_inv,
511 col_perm,
512 col_perm_inv,
513 par.reborrow(),
514 stack,
515 params,
516 );
517 linalg::copy_across_layouts(mat, work.rb(), par);
518 return ret;
519 }
520 let size = Ord::min(M, N);
521 let (mut row_transpositions, stack) = stack.make_with(size, |_| 0);
522 let row_transpositions = row_transpositions.as_mut();
523 let (mut col_transpositions, _) = stack.make_with(size, |_| 0);
524 let col_transpositions = col_transpositions.as_mut();
525 let n_transpositions = if mat.row_stride().abs() < mat.col_stride().abs() {
526 lu_in_place_unblocked(
527 mat,
528 row_transpositions,
529 col_transpositions,
530 par,
531 false,
532 params,
533 )
534 } else {
535 lu_in_place_unblocked(
536 mat.transpose_mut(),
537 col_transpositions,
538 row_transpositions,
539 par,
540 true,
541 params,
542 )
543 };
544 for i in 0..M {
545 row_perm[i] = I::truncate(i);
546 }
547 for (i, t) in row_transpositions.iter().copied().enumerate() {
548 row_perm.as_mut().swap(i, t);
549 }
550 for i in 0..M {
551 row_perm_inv[row_perm[i].zx()] = I::truncate(i);
552 }
553 for j in 0..N {
554 col_perm[j] = I::truncate(j);
555 }
556 for (i, t) in col_transpositions.iter().copied().enumerate() {
557 col_perm.as_mut().swap(i, t);
558 }
559 for j in 0..N {
560 col_perm_inv[col_perm[j].zx()] = I::truncate(j);
561 }
562 unsafe {
563 (
564 FullPivLuInfo {
565 transposition_count: n_transpositions,
566 },
567 PermRef::new_unchecked(row_perm, row_perm_inv, M),
568 PermRef::new_unchecked(col_perm, col_perm_inv, N),
569 )
570 }
571}
572#[cfg(test)]
573mod tests {
574 use super::*;
575 use crate::stats::prelude::*;
576 use crate::utils::approx::*;
577 use crate::{Mat, assert, c64};
578 use dyn_stack::MemBuffer;
579 #[test]
580 fn test_flu() {
581 let rng = &mut StdRng::seed_from_u64(0);
582 for mut par in [Par::Seq, Par::rayon(8)] {
583 let par_shape = par.shape();
584 for m in [8, 16, 24, 32, 128, 255, 256, 257] {
585 let n = 8;
586 let approx_eq = CwiseMat(ApproxEq {
587 abs_tol: 1e-10,
588 rel_tol: 1e-10,
589 });
590 let A = CwiseMatDistribution {
591 nrows: m,
592 ncols: n,
593 dist: ComplexDistribution::new(
594 StandardNormal,
595 StandardNormal,
596 ),
597 }
598 .rand::<Mat<c64>>(rng);
599 let A = A.as_ref();
600 let mut LU = A.cloned();
601 let row_perm = &mut *vec![0usize; m];
602 let row_perm_inv = &mut *vec![0usize; m];
603 let col_perm = &mut *vec![0usize; n];
604 let col_perm_inv = &mut *vec![0usize; n];
605 let (_, p, q) = lu_in_place(
606 LU.as_mut(),
607 row_perm,
608 row_perm_inv,
609 col_perm,
610 col_perm_inv,
611 par.reborrow(),
612 MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
613 usize,
614 c64,
615 >(
616 n, n, par_shape, default()
617 ))),
618 default(),
619 );
620 let mut L = LU.as_ref().cloned();
621 let mut U = LU.as_ref().cloned();
622 for j in 0..n {
623 for i in 0..j {
624 L[(i, j)] = c64::ZERO;
625 }
626 L[(j, j)] = c64::ONE;
627 }
628 for j in 0..n {
629 for i in j + 1..m {
630 U[(i, j)] = c64::ZERO;
631 }
632 }
633 let L = L.as_ref();
634 let U = U.as_ref();
635 let U = U.subrows(0, n);
636 assert!(p.inverse() * L * U * q ~ A);
637 }
638 for n in [16, 24, 32, 128, 255, 256, 257] {
639 let approx_eq = CwiseMat(ApproxEq {
640 abs_tol: 1e-10,
641 rel_tol: 1e-10,
642 });
643 let A = CwiseMatDistribution {
644 nrows: n,
645 ncols: n,
646 dist: StandardNormal,
647 }
648 .rand::<Mat<f64>>(rng);
649 let A = A.as_ref();
650 let mut LU = A.cloned();
651 let row_perm = &mut *vec![0usize; n];
652 let row_perm_inv = &mut *vec![0usize; n];
653 let col_perm = &mut *vec![0usize; n];
654 let col_perm_inv = &mut *vec![0usize; n];
655 let (_, p, q) = lu_in_place(
656 LU.as_mut(),
657 row_perm,
658 row_perm_inv,
659 col_perm,
660 col_perm_inv,
661 par.reborrow(),
662 MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
663 usize,
664 f64,
665 >(
666 n, n, par_shape, default()
667 ))),
668 default(),
669 );
670 let mut L = LU.as_ref().cloned();
671 let mut U = LU.as_ref().cloned();
672 for j in 0..n {
673 for i in 0..j {
674 L[(i, j)] = 0.0;
675 }
676 L[(j, j)] = 1.0;
677 }
678 for j in 0..n {
679 for i in j + 1..n {
680 U[(i, j)] = 0.0;
681 }
682 }
683 let L = L.as_ref();
684 let U = U.as_ref();
685 assert!(p.inverse() * L * U * q ~ A);
686 }
687 }
688 }
689 #[test]
690 fn test_gh238() {
691 #[rustfmt::skip]
692 let A = [
693 [
694 -2.0,
695 0.0,
696 1.0,
697 0.0,
698 0.0,
699 0.0,
700 0.0,
701 0.0,
702 0.0,
703 0.0,
704 0.0,
705 0.0,
706 0.0,
707 0.0,
708 0.0,
709 0.0,
710 0.0,
711 0.0,
712 0.0,
713 0.0,
714 0.0,
715 0.0,
716 0.0,
717 0.0,
718 0.0,
719 0.0,
720 0.0,
721 0.0,
722 ],
723 [
724 0.0,
725 -2.0,
726 0.0,
727 1.0,
728 0.0,
729 0.0,
730 0.0,
731 0.0,
732 0.0,
733 0.0,
734 0.0,
735 0.0,
736 0.0,
737 0.0,
738 0.0,
739 0.0,
740 0.0,
741 0.0,
742 0.0,
743 0.0,
744 0.0,
745 0.0,
746 0.0,
747 0.0,
748 0.0,
749 0.0,
750 0.0,
751 0.0,
752 ],
753 [
754 0.0,
755 0.0,
756 0.0,
757 0.0,
758 -2.0,
759 0.0,
760 1.0,
761 0.0,
762 0.0,
763 0.0,
764 0.0,
765 0.0,
766 0.0,
767 0.0,
768 0.0,
769 0.0,
770 0.0,
771 0.0,
772 0.0,
773 0.0,
774 0.0,
775 0.0,
776 0.0,
777 0.0,
778 0.0,
779 0.0,
780 0.0,
781 0.0,
782 ],
783 [
784 0.0,
785 0.0,
786 0.0,
787 0.0,
788 0.0,
789 -2.0,
790 0.0,
791 1.0,
792 0.0,
793 0.0,
794 0.0,
795 0.0,
796 0.0,
797 0.0,
798 0.0,
799 0.0,
800 0.0,
801 0.0,
802 0.0,
803 0.0,
804 0.0,
805 0.0,
806 0.0,
807 0.0,
808 0.0,
809 0.0,
810 0.0,
811 0.0,
812 ],
813 [
814 0.0,
815 0.0,
816 0.0,
817 0.0,
818 0.0,
819 0.0,
820 0.0,
821 0.0,
822 -2.0,
823 0.0,
824 1.0,
825 0.0,
826 0.0,
827 0.0,
828 0.0,
829 0.0,
830 0.0,
831 0.0,
832 0.0,
833 0.0,
834 0.0,
835 0.0,
836 0.0,
837 0.0,
838 0.0,
839 0.0,
840 0.0,
841 0.0,
842 ],
843 [
844 0.0,
845 0.0,
846 0.0,
847 0.0,
848 0.0,
849 0.0,
850 0.0,
851 0.0,
852 0.0,
853 -2.0,
854 0.0,
855 1.0,
856 0.0,
857 0.0,
858 0.0,
859 0.0,
860 0.0,
861 0.0,
862 0.0,
863 0.0,
864 0.0,
865 0.0,
866 0.0,
867 0.0,
868 0.0,
869 0.0,
870 0.0,
871 0.0,
872 ],
873 [
874 0.0,
875 0.0,
876 0.0,
877 0.0,
878 0.0,
879 0.0,
880 0.0,
881 0.0,
882 0.0,
883 0.0,
884 0.0,
885 0.0,
886 -2.0,
887 0.0,
888 1.0,
889 0.0,
890 0.0,
891 0.0,
892 0.0,
893 0.0,
894 0.0,
895 0.0,
896 0.0,
897 0.0,
898 0.0,
899 0.0,
900 0.0,
901 0.0,
902 ],
903 [
904 0.0,
905 0.0,
906 0.0,
907 0.0,
908 0.0,
909 0.0,
910 0.0,
911 0.0,
912 0.0,
913 0.0,
914 0.0,
915 0.0,
916 0.0,
917 -2.0,
918 0.0,
919 1.0,
920 0.0,
921 0.0,
922 0.0,
923 0.0,
924 0.0,
925 0.0,
926 0.0,
927 0.0,
928 0.0,
929 0.0,
930 0.0,
931 0.0,
932 ],
933 [
934 0.0,
935 0.0,
936 0.0,
937 0.0,
938 0.0,
939 0.0,
940 0.0,
941 0.0,
942 0.0,
943 0.0,
944 0.0,
945 0.0,
946 0.0,
947 0.0,
948 0.0,
949 0.0,
950 -2.0,
951 0.0,
952 1.0,
953 0.0,
954 0.0,
955 0.0,
956 0.0,
957 0.0,
958 0.0,
959 0.0,
960 0.0,
961 0.0,
962 ],
963 [
964 0.0,
965 0.0,
966 0.0,
967 0.0,
968 0.0,
969 0.0,
970 0.0,
971 0.0,
972 0.0,
973 0.0,
974 0.0,
975 0.0,
976 0.0,
977 0.0,
978 0.0,
979 0.0,
980 0.0,
981 -2.0,
982 0.0,
983 1.0,
984 0.0,
985 0.0,
986 0.0,
987 0.0,
988 0.0,
989 0.0,
990 0.0,
991 0.0,
992 ],
993 [
994 0.0,
995 0.0,
996 0.0,
997 0.0,
998 0.0,
999 0.0,
1000 0.0,
1001 0.0,
1002 0.0,
1003 0.0,
1004 0.0,
1005 0.0,
1006 0.0,
1007 0.0,
1008 0.0,
1009 0.0,
1010 0.0,
1011 0.0,
1012 0.0,
1013 0.0,
1014 -2.0,
1015 0.0,
1016 1.0,
1017 0.0,
1018 0.0,
1019 0.0,
1020 0.0,
1021 0.0,
1022 ],
1023 [
1024 0.0,
1025 0.0,
1026 0.0,
1027 0.0,
1028 0.0,
1029 0.0,
1030 0.0,
1031 0.0,
1032 0.0,
1033 0.0,
1034 0.0,
1035 0.0,
1036 0.0,
1037 0.0,
1038 0.0,
1039 0.0,
1040 0.0,
1041 0.0,
1042 0.0,
1043 0.0,
1044 0.0,
1045 -2.0,
1046 0.0,
1047 1.0,
1048 0.0,
1049 0.0,
1050 0.0,
1051 0.0,
1052 ],
1053 [
1054 0.0,
1055 0.0,
1056 0.0,
1057 0.0,
1058 0.0,
1059 0.0,
1060 0.0,
1061 0.0,
1062 0.0,
1063 0.0,
1064 0.0,
1065 0.0,
1066 0.0,
1067 0.0,
1068 0.0,
1069 0.0,
1070 0.0,
1071 0.0,
1072 0.0,
1073 0.0,
1074 0.0,
1075 0.0,
1076 0.0,
1077 0.0,
1078 -2.0,
1079 0.0,
1080 1.0,
1081 0.0,
1082 ],
1083 [
1084 0.0,
1085 0.0,
1086 0.0,
1087 0.0,
1088 0.0,
1089 0.0,
1090 0.0,
1091 0.0,
1092 0.0,
1093 0.0,
1094 0.0,
1095 0.0,
1096 0.0,
1097 0.0,
1098 0.0,
1099 0.0,
1100 0.0,
1101 0.0,
1102 0.0,
1103 0.0,
1104 0.0,
1105 0.0,
1106 0.0,
1107 0.0,
1108 0.0,
1109 -2.0,
1110 0.0,
1111 1.0,
1112 ],
1113 [
1114 0.0,
1115 0.0,
1116 0.0,
1117 0.0,
1118 0.0,
1119 0.0,
1120 0.0,
1121 0.0,
1122 0.0,
1123 0.0,
1124 0.0,
1125 0.0,
1126 0.0,
1127 0.0,
1128 0.0,
1129 0.0,
1130 0.0,
1131 0.0,
1132 0.0,
1133 0.0,
1134 0.0,
1135 0.0,
1136 0.0,
1137 0.0,
1138 0.0,
1139 0.0,
1140 0.0,
1141 1.0,
1142 ],
1143 [
1144 0.0,
1145 0.0,
1146 0.0,
1147 1.0,
1148 0.0,
1149 0.0,
1150 1.0,
1151 0.0,
1152 0.0,
1153 0.0,
1154 0.0,
1155 0.0,
1156 0.0,
1157 0.0,
1158 0.0,
1159 0.0,
1160 0.0,
1161 0.0,
1162 0.0,
1163 0.0,
1164 0.0,
1165 0.0,
1166 0.0,
1167 0.0,
1168 0.0,
1169 0.0,
1170 0.0,
1171 0.0,
1172 ],
1173 [
1174 0.0,
1175 -1.0,
1176 0.0,
1177 0.0,
1178 1.0,
1179 0.0,
1180 0.0,
1181 0.0,
1182 0.0,
1183 0.0,
1184 0.0,
1185 0.0,
1186 0.0,
1187 0.0,
1188 0.0,
1189 0.0,
1190 0.0,
1191 0.0,
1192 0.0,
1193 0.0,
1194 0.0,
1195 0.0,
1196 0.0,
1197 0.0,
1198 0.0,
1199 0.0,
1200 0.0,
1201 0.0,
1202 ],
1203 [
1204 0.0,
1205 0.0,
1206 0.0,
1207 0.0,
1208 0.0,
1209 0.0,
1210 0.0,
1211 0.0,
1212 0.0,
1213 0.0,
1214 0.0,
1215 0.0,
1216 0.0,
1217 0.0,
1218 0.0,
1219 1.0,
1220 0.0,
1221 0.0,
1222 1.0,
1223 0.0,
1224 0.0,
1225 0.0,
1226 0.0,
1227 0.0,
1228 0.0,
1229 0.0,
1230 0.0,
1231 0.0,
1232 ],
1233 [
1234 0.0,
1235 0.0,
1236 0.0,
1237 0.0,
1238 0.0,
1239 0.0,
1240 0.0,
1241 0.0,
1242 0.0,
1243 0.0,
1244 0.0,
1245 0.0,
1246 0.0,
1247 -1.0,
1248 0.0,
1249 0.0,
1250 1.0,
1251 0.0,
1252 0.0,
1253 0.0,
1254 0.0,
1255 0.0,
1256 0.0,
1257 0.0,
1258 0.0,
1259 0.0,
1260 0.0,
1261 0.0,
1262 ],
1263 [
1264 0.0,
1265 0.0,
1266 0.0,
1267 0.0,
1268 0.0,
1269 0.0,
1270 0.0,
1271 1.0,
1272 0.0,
1273 0.0,
1274 1.0,
1275 0.0,
1276 0.0,
1277 0.0,
1278 0.0,
1279 0.0,
1280 0.0,
1281 0.0,
1282 0.0,
1283 0.0,
1284 0.0,
1285 0.0,
1286 0.0,
1287 0.0,
1288 0.0,
1289 0.0,
1290 0.0,
1291 0.0,
1292 ],
1293 [
1294 0.0,
1295 0.0,
1296 0.0,
1297 0.0,
1298 0.0,
1299 1.0,
1300 0.0,
1301 0.0,
1302 -1.0,
1303 0.0,
1304 0.0,
1305 0.0,
1306 0.0,
1307 0.0,
1308 0.0,
1309 0.0,
1310 0.0,
1311 0.0,
1312 0.0,
1313 0.0,
1314 0.0,
1315 0.0,
1316 0.0,
1317 0.0,
1318 0.0,
1319 0.0,
1320 0.0,
1321 0.0,
1322 ],
1323 [
1324 0.0,
1325 0.0,
1326 0.0,
1327 0.0,
1328 0.0,
1329 0.0,
1330 0.0,
1331 0.0,
1332 0.0,
1333 0.0,
1334 0.0,
1335 1.0,
1336 0.0,
1337 0.0,
1338 1.0,
1339 0.0,
1340 0.0,
1341 0.0,
1342 0.0,
1343 0.0,
1344 0.0,
1345 0.0,
1346 0.0,
1347 0.0,
1348 0.0,
1349 0.0,
1350 0.0,
1351 0.0,
1352 ],
1353 [
1354 0.0,
1355 0.0,
1356 0.0,
1357 0.0,
1358 0.0,
1359 0.0,
1360 0.0,
1361 0.0,
1362 0.0,
1363 1.0,
1364 0.0,
1365 0.0,
1366 -1.0,
1367 0.0,
1368 0.0,
1369 0.0,
1370 0.0,
1371 0.0,
1372 0.0,
1373 0.0,
1374 0.0,
1375 0.0,
1376 0.0,
1377 0.0,
1378 0.0,
1379 0.0,
1380 0.0,
1381 0.0,
1382 ],
1383 [
1384 0.0,
1385 0.0,
1386 0.0,
1387 0.0,
1388 0.0,
1389 0.0,
1390 0.0,
1391 0.0,
1392 0.0,
1393 0.0,
1394 0.0,
1395 0.0,
1396 0.0,
1397 0.0,
1398 0.0,
1399 0.0,
1400 0.0,
1401 0.0,
1402 0.0,
1403 0.0,
1404 0.0,
1405 0.0,
1406 0.0,
1407 1.0,
1408 0.0,
1409 0.0,
1410 1.0,
1411 0.0,
1412 ],
1413 [
1414 0.0,
1415 0.0,
1416 0.0,
1417 0.0,
1418 0.0,
1419 0.0,
1420 0.0,
1421 0.0,
1422 0.0,
1423 0.0,
1424 0.0,
1425 0.0,
1426 0.0,
1427 0.0,
1428 0.0,
1429 0.0,
1430 0.0,
1431 0.0,
1432 0.0,
1433 0.0,
1434 0.0,
1435 -1.0,
1436 0.0,
1437 0.0,
1438 1.0,
1439 0.0,
1440 0.0,
1441 0.0,
1442 ],
1443 [
1444 0.0,
1445 0.0,
1446 1.0,
1447 0.0,
1448 0.0,
1449 0.0,
1450 0.0,
1451 0.0,
1452 0.0,
1453 0.0,
1454 0.0,
1455 0.0,
1456 0.0,
1457 0.0,
1458 0.0,
1459 0.0,
1460 0.0,
1461 0.0,
1462 0.0,
1463 0.0,
1464 0.0,
1465 0.0,
1466 0.0,
1467 0.0,
1468 0.0,
1469 0.0,
1470 0.0,
1471 0.0,
1472 ],
1473 [
1474 0.0,
1475 0.0,
1476 0.0,
1477 0.0,
1478 0.0,
1479 0.0,
1480 0.0,
1481 0.0,
1482 0.0,
1483 0.0,
1484 0.0,
1485 0.0,
1486 0.0,
1487 0.0,
1488 0.0,
1489 0.0,
1490 0.0,
1491 0.0,
1492 0.0,
1493 1.0,
1494 0.0,
1495 0.0,
1496 1.0,
1497 0.0,
1498 0.0,
1499 0.0,
1500 0.0,
1501 0.0,
1502 ],
1503 [
1504 0.0,
1505 0.0,
1506 0.0,
1507 0.0,
1508 0.0,
1509 0.0,
1510 0.0,
1511 0.0,
1512 0.0,
1513 0.0,
1514 0.0,
1515 0.0,
1516 0.0,
1517 0.0,
1518 0.0,
1519 0.0,
1520 0.0,
1521 1.0,
1522 0.0,
1523 0.0,
1524 -1.0,
1525 0.0,
1526 0.0,
1527 0.0,
1528 0.0,
1529 0.0,
1530 0.0,
1531 0.0,
1532 ],
1533 ];
1534 let A = MatRef::from_row_major_array(&A).to_owned();
1535 let b = [
1536 -0.0, -0.0, -2.0, -2.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0,
1537 -0.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1538 0.0, 0.0, 0.0, 0.0,
1539 ];
1540 let b = ColRef::from_slice(&b);
1541 let lu = A.full_piv_lu(Par::Seq);
1542 assert!(lu.solve(b, Par::Seq).is_all_finite());
1543 }
1544}