1use crate::internal_prelude::*;
2use crate::perm::{swap_cols_idx, swap_rows_idx};
3use crate::utils::thread::par_split_indices;
4use faer_traits::{Real, RealReg};
5use linalg::matmul::matmul;
6use pulp::Simd;
7
8#[inline(always)]
9fn best_value<T: ComplexField, S: Simd>(
10 simd: &SimdCtx<T, S>,
11 best_value: RealReg<T::SimdVec<S>>,
12 best_indices: T::SimdIndex<S>,
13 value: T::SimdVec<S>,
14 indices: T::SimdIndex<S>,
15) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
16 let value = simd.abs1(value);
17 let is_better = (**simd).gt(value, best_value);
18 (
19 RealReg(simd.select(is_better, value.0, best_value.0)),
20 simd.iselect(is_better, indices, best_indices),
21 )
22}
23#[inline(always)]
24fn best_score<T: ComplexField, S: Simd>(
25 simd: &SimdCtx<T, S>,
26 best_score: RealReg<T::SimdVec<S>>,
27 best_indices: T::SimdIndex<S>,
28 score: RealReg<T::SimdVec<S>>,
29 indices: T::SimdIndex<S>,
30) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
31 let is_better = (**simd).gt(score, best_score);
32 (
33 RealReg(simd.select(is_better, score.0, best_score.0)),
34 simd.iselect(is_better, indices, best_indices),
35 )
36}
37#[inline(always)]
38fn best_score_2d<T: ComplexField, S: Simd>(
39 simd: &SimdCtx<T, S>,
40 best_score: RealReg<T::SimdVec<S>>,
41 best_row: T::SimdIndex<S>,
42 best_col: T::SimdIndex<S>,
43 score: RealReg<T::SimdVec<S>>,
44 row: T::SimdIndex<S>,
45 col: T::SimdIndex<S>,
46) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>, T::SimdIndex<S>) {
47 let is_better = (**simd).gt(score, best_score);
48 (
49 RealReg(simd.select(is_better, score.0, best_score.0)),
50 simd.iselect(is_better, row, best_row),
51 simd.iselect(is_better, col, best_col),
52 )
53}
54#[inline(always)]
55fn reduce_2d<T: ComplexField, S: Simd>(
56 simd: &SimdCtx<T, S>,
57 best_values: RealReg<T::SimdVec<S>>,
58 best_row: T::SimdIndex<S>,
59 best_col: T::SimdIndex<S>,
60) -> (usize, usize, Real<T>) {
61 let best_val = simd.reduce_max_real(best_values);
62 let best_val_splat = simd.splat_real(&best_val);
63 let is_best = (**simd).ge(best_values, best_val_splat);
64 let idx = simd.first_true_mask(is_best);
65 let best_row = bytemuck::cast_slice::<T::SimdIndex<S>, T::Index>(
66 core::slice::from_ref(&best_row),
67 )[idx];
68 let best_col = bytemuck::cast_slice::<T::SimdIndex<S>, T::Index>(
69 core::slice::from_ref(&best_col),
70 )[idx];
71 (best_row.zx(), best_col.zx(), best_val)
72}
73#[inline(always)]
74fn best_in_col_simd<'M, T: ComplexField, S: Simd>(
75 simd: SimdCtx<'M, T, S>,
76 data: ColRef<'_, T, Dim<'M>, ContiguousFwd>,
77) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
78 let iota = T::simd_iota(&simd.0);
79 let lane_count =
80 core::mem::size_of::<T::SimdVec<S>>() / core::mem::size_of::<T>();
81 let inc = simd.isplat(T::Index::truncate(lane_count));
82
83 let indices = simd.batch_indices::<4>();
84 let mut best_val = [RealReg(simd.zero()); 4];
85 let mut best_idx = [simd.izero(); 4];
86
87 let mut idx = simd.iadd(
88 iota,
89 simd.isplat(T::Index::truncate(simd.offset().wrapping_neg())),
90 );
91 simd_iter!(for (IDX, i) in [indices; 4] {
92 (best_val[IDX], best_idx[IDX]) = best_value(
93 &simd,
94 best_val[IDX],
95 best_idx[IDX],
96 simd.read(data, i),
97 idx,
98 );
99 idx = simd.iadd(idx, inc);
100 });
101 let _ = idx;
102
103 let (best_val0, best_idx0) =
104 best_score(&simd, best_val[0], best_idx[0], best_val[1], best_idx[1]);
105 let (best_val2, best_idx2) =
106 best_score(&simd, best_val[2], best_idx[2], best_val[3], best_idx[3]);
107 best_score(&simd, best_val0, best_idx0, best_val2, best_idx2)
108}
109#[inline(always)]
110fn update_and_best_in_col_simd<'M, T: ComplexField, S: Simd>(
111 simd: SimdCtx<'M, T, S>,
112 data: ColMut<'_, T, Dim<'M>, ContiguousFwd>,
113 lhs: ColRef<'_, T, Dim<'M>, ContiguousFwd>,
114 rhs: T,
115) -> (RealReg<T::SimdVec<S>>, T::SimdIndex<S>) {
116 let mut data = data;
117 let iota = T::simd_iota(&simd.0);
118 let lane_count =
119 core::mem::size_of::<T::SimdVec<S>>() / core::mem::size_of::<T>();
120 let inc = simd.isplat(T::Index::truncate(lane_count));
121 let mut best_val = [RealReg(simd.zero()); 3];
122 let mut best_idx = [simd.izero(); 3];
123
124 let mut idx = simd.iadd(
125 iota,
126 simd.isplat(T::Index::truncate(simd.offset().wrapping_neg())),
127 );
128
129 let indices = simd.batch_indices::<3>();
130 let rhs = simd.splat(&-rhs);
131 simd_iter!(for (IDX, i) in [indices; 3] {
132 let mut x0 = simd.read(data.rb(), i);
133 let l0 = simd.read(lhs, i);
134 x0 = simd.mul_add(l0, rhs, x0);
135 (best_val[IDX], best_idx[IDX]) =
136 best_value(&simd, best_val[IDX], best_idx[IDX], x0, idx);
137 simd.write(data.rb_mut(), i, x0);
138 idx = simd.iadd(idx, inc);
139 });
140 let _ = idx;
141
142 let (best_val0, best_idx0) =
143 best_score(&simd, best_val[0], best_idx[0], best_val[1], best_idx[1]);
144 let (best_val2, best_idx2) = (best_val[2], best_idx[2]);
145 best_score(&simd, best_val0, best_idx0, best_val2, best_idx2)
146}
147#[inline(always)]
148fn best_in_mat_simd<T: ComplexField>(
149 data: MatRef<'_, T, usize, usize, ContiguousFwd>,
150) -> (usize, usize, Real<T>) {
151 struct Impl<'a, 'M, 'N, T: ComplexField> {
152 data: MatRef<'a, T, Dim<'M>, Dim<'N>, ContiguousFwd>,
153 }
154 impl<'a, 'M, 'N, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'N, T> {
155 type Output = (usize, usize, Real<T>);
156
157 #[inline(always)]
158 fn with_simd<S: Simd>(self, simd: S) -> Self::Output {
159 let Self { data } = self;
160 let M = data.nrows();
161 let N = data.ncols();
162 let simd = SimdCtx::<'_, T, S>::new(T::simd_ctx(simd), M);
163 let mut best_row = simd.isplat(T::Index::truncate(0));
164 let mut best_col = simd.isplat(T::Index::truncate(0));
165 let mut best_val = simd.splat_real(&zero());
166 for j in N.indices() {
167 let col = data.col(j);
168 let (best_val_j, best_row_j) = best_in_col_simd(simd, col);
169 (best_val, best_row, best_col) = best_score_2d(
170 &simd,
171 best_val,
172 best_row,
173 best_col,
174 best_val_j,
175 best_row_j,
176 simd.isplat(T::Index::truncate(*j)),
177 );
178 }
179 reduce_2d(&simd, best_val, best_row, best_col)
180 }
181 }
182 with_dim!(M, data.nrows());
183 with_dim!(N, data.ncols());
184 dispatch!(
185 Impl {
186 data: data.as_shape(M, N)
187 },
188 Impl,
189 T
190 )
191}
192#[inline(always)]
193fn update_and_best_in_mat_simd<T: ComplexField>(
194 data: MatMut<'_, T, usize, usize, ContiguousFwd>,
195 lhs: ColRef<'_, T, usize, ContiguousFwd>,
196 rhs: RowRef<'_, T, usize>,
197 align: usize,
198) -> (usize, usize, Real<T>) {
199 struct Impl<'a, 'M, 'N, T: ComplexField> {
200 data: MatMut<'a, T, Dim<'M>, Dim<'N>, ContiguousFwd>,
201 lhs: ColRef<'a, T, Dim<'M>, ContiguousFwd>,
202 rhs: RowRef<'a, T, Dim<'N>>,
203 align: usize,
204 }
205 impl<'a, 'M, 'N, T: ComplexField> pulp::WithSimd for Impl<'a, 'M, 'N, T> {
206 type Output = (usize, usize, Real<T>);
207
208 #[inline(always)]
209 fn with_simd<S: Simd>(self, simd: S) -> Self::Output {
210 let Self {
211 data,
212 lhs,
213 rhs,
214 align,
215 } = self;
216 let M = data.nrows();
217 let N = data.ncols();
218 let simd =
219 SimdCtx::<'_, T, S>::new_align(T::simd_ctx(simd), M, align);
220 let mut best_row = simd.isplat(T::Index::truncate(0));
221 let mut best_col = simd.isplat(T::Index::truncate(0));
222 let mut best_val = simd.splat_real(&zero());
223 let mut data = data;
224 for j in N.indices() {
225 let data = data.rb_mut().col_mut(j);
226 let rhs = rhs[j].copy();
227 let (best_val_j, best_row_j) =
228 update_and_best_in_col_simd(simd, data, lhs, rhs);
229 (best_val, best_row, best_col) = best_score_2d(
230 &simd,
231 best_val,
232 best_row,
233 best_col,
234 best_val_j,
235 best_row_j,
236 simd.isplat(T::Index::truncate(*j)),
237 );
238 }
239 reduce_2d(&simd, best_val, best_row, best_col)
240 }
241 }
242 with_dim!(M, data.nrows());
243 with_dim!(N, data.ncols());
244 dispatch!(
245 Impl {
246 data: data.as_shape_mut(M, N),
247 lhs: lhs.as_row_shape(M),
248 rhs: rhs.as_col_shape(N),
249 align,
250 },
251 Impl,
252 T
253 )
254}
255fn best_in_matrix_fallback<T: ComplexField>(
256 data: MatRef<'_, T>,
257) -> (usize, usize, Real<T>) {
258 let mut max = zero();
259 let mut row = 0;
260 let mut col = 0;
261 let (m, n) = data.shape();
262 for j in 0..n {
263 for i in 0..m {
264 let abs = data[(i, j)].abs1();
265 if abs > max {
266 row = i;
267 col = j;
268 max = abs;
269 }
270 }
271 }
272 (row, col, max)
273}
274fn best_in_matrix<T: ComplexField>(
275 data: MatRef<'_, T>,
276) -> (usize, usize, Real<T>) {
277 if const { T::SIMD_CAPABILITIES.is_simd() } {
278 if let Some(dst) = data.try_as_col_major() {
279 best_in_mat_simd(dst)
280 } else {
281 best_in_matrix_fallback(data)
282 }
283 } else {
284 best_in_matrix_fallback(data)
285 }
286}
287fn rank_one_update_and_best_in_matrix<T: ComplexField>(
288 mut dst: MatMut<'_, T>,
289 lhs: ColRef<'_, T>,
290 rhs: RowRef<'_, T>,
291 align: usize,
292) -> (usize, usize, Real<T>) {
293 if const { T::SIMD_CAPABILITIES.is_simd() } {
294 if let (Some(dst), Some(lhs)) =
295 (dst.rb_mut().try_as_col_major_mut(), lhs.try_as_col_major())
296 {
297 update_and_best_in_mat_simd(dst, lhs, rhs, align)
298 } else {
299 matmul(
300 dst.rb_mut(),
301 Accum::Add,
302 lhs.as_mat(),
303 rhs.as_mat(),
304 -one::<T>(),
305 Par::Seq,
306 );
307 best_in_matrix(dst.rb())
308 }
309 } else {
310 matmul(
311 dst.rb_mut(),
312 Accum::Add,
313 lhs.as_mat(),
314 rhs.as_mat(),
315 -one::<T>(),
316 Par::Seq,
317 );
318 best_in_matrix(dst.rb())
319 }
320}
321fn lu_in_place_unblocked<T: ComplexField>(
322 A: MatMut<'_, T>,
323 row_trans: &mut [usize],
324 col_trans: &mut [usize],
325 par: Par,
326 transpose: bool,
327 params: Spec<FullPivLuParams, T>,
328) -> usize {
329 let params = params.config;
330 let mut n_trans = 0;
331 let (m, n) = A.shape();
332 if m == 0 || n == 0 {
333 return 0;
334 }
335 let mut par = par;
336 let mut A = A;
337 let (mut max_row, mut max_col, mut max_score) = best_in_matrix(A.rb());
338 for k in 0..Ord::min(m, n) {
339 if max_score < min_positive() {
340 for (i, (row, col)) in
341 core::iter::zip(&mut row_trans[k..], &mut col_trans[k..])
342 .enumerate()
343 {
344 *row = i + k;
345 *col = i + k;
346 }
347 break;
348 }
349 row_trans[k] = max_row;
350 col_trans[k] = max_col;
351 if max_row != k {
352 swap_rows_idx(A.rb_mut(), k, max_row);
353 n_trans += 1;
354 }
355 if max_col != k {
356 swap_cols_idx(A.rb_mut(), k, max_col);
357 n_trans += 1;
358 }
359 let ref inv = A[(k, k)].recip();
360 if transpose {
361 for j in k + 1..n {
362 A[(k, j)] *= inv;
363 }
364 } else {
365 for i in k + 1..m {
366 A[(i, k)] *= inv;
367 }
368 }
369 if k + 1 == Ord::min(m, n) {
370 break;
371 }
372 if (m - k - 1) * (n - k - 1) < params.par_threshold {
373 par = Par::Seq;
374 }
375 let (_, A01, A10, mut A11) = A.rb_mut().split_at_mut(k + 1, k + 1);
376 let lhs = A10.col(k);
377 let rhs = A01.row(k);
378 match par {
379 Par::Seq => {
380 (max_row, max_col, max_score) =
381 rank_one_update_and_best_in_matrix(
382 A11.rb_mut(),
383 lhs,
384 rhs,
385 simd_align(k + 1),
386 );
387 },
388 #[cfg(feature = "rayon")]
389 Par::Rayon(nthreads) => {
390 use rayon::prelude::*;
391 let nthreads = nthreads.get();
392 let mut best = core::iter::repeat_with(|| (0, 0, zero()))
393 .take(nthreads)
394 .collect::<alloc::vec::Vec<_>>();
395 let full_cols = A11.ncols();
396 spindle::for_each(
397 nthreads,
398 best.par_iter_mut()
399 .zip_eq(A11.rb_mut().par_col_partition_mut(nthreads))
400 .zip_eq(rhs.par_partition(nthreads))
401 .enumerate(),
402 |(idx, (((max_row, max_col, max_score), A11), rhs))| {
403 (*max_row, *max_col, *max_score) = {
404 let (a, mut b, c) =
405 rank_one_update_and_best_in_matrix(
406 A11,
407 lhs,
408 rhs,
409 simd_align(k + 1),
410 );
411 b += par_split_indices(full_cols, idx, nthreads).0;
412 (a, b, c)
413 };
414 },
415 );
416 max_row = 0;
417 max_col = 0;
418 max_score = zero();
419 for (row, col, val) in best {
420 if val > max_score {
421 max_row = row;
422 max_col = col;
423 max_score = val;
424 }
425 }
426 },
427 }
428 max_row += k + 1;
429 max_col += k + 1;
430 }
431 n_trans
432}
433#[derive(Copy, Clone, Debug)]
435pub struct FullPivLuParams {
436 pub par_threshold: usize,
438 #[doc(hidden)]
439 pub non_exhaustive: NonExhaustive,
440}
441impl<T: ComplexField> Auto<T> for FullPivLuParams {
442 #[inline]
443 fn auto() -> Self {
444 Self {
445 par_threshold: 256 * 512,
446 non_exhaustive: NonExhaustive(()),
447 }
448 }
449}
450#[inline]
451pub fn lu_in_place_scratch<I: Index, T: ComplexField>(
452 nrows: usize,
453 ncols: usize,
454 par: Par,
455 params: Spec<FullPivLuParams, T>,
456) -> StackReq {
457 _ = par;
458 _ = params;
459 let size = Ord::min(nrows, ncols);
460 StackReq::new::<usize>(size).array(2)
461}
462#[derive(Copy, Clone, Debug)]
463pub struct FullPivLuInfo {
464 pub transposition_count: usize,
465}
466pub fn lu_in_place<'out, I: Index, T: ComplexField>(
467 mat: MatMut<'_, T>,
468 row_perm: &'out mut [I],
469 row_perm_inv: &'out mut [I],
470 col_perm: &'out mut [I],
471 col_perm_inv: &'out mut [I],
472 par: Par,
473 stack: &mut MemStack,
474 params: Spec<FullPivLuParams, T>,
475) -> (FullPivLuInfo, PermRef<'out, I>, PermRef<'out, I>) {
476 #[cfg(feature = "perf-warn")]
477 if (mat.col_stride().unsigned_abs() == 1
478 || mat.row_stride().unsigned_abs() != 1)
479 && crate::__perf_warn!(LU_WARN)
480 {
481 log::warn!(
482 target : "faer_perf",
483 "LU with full pivoting prefers column-major or row-major matrix. Found matrix with generic strides."
484 );
485 }
486 let (M, N) = mat.shape();
487 let size = Ord::min(M, N);
488 let (mut row_transpositions, stack) = stack.make_with(size, |_| 0);
489 let row_transpositions = row_transpositions.as_mut();
490 let (mut col_transpositions, _) = stack.make_with(size, |_| 0);
491 let col_transpositions = col_transpositions.as_mut();
492 let n_transpositions = if mat.row_stride().abs() < mat.col_stride().abs() {
493 lu_in_place_unblocked(
494 mat,
495 row_transpositions,
496 col_transpositions,
497 par,
498 false,
499 params,
500 )
501 } else {
502 lu_in_place_unblocked(
503 mat.transpose_mut(),
504 col_transpositions,
505 row_transpositions,
506 par,
507 true,
508 params,
509 )
510 };
511 for i in 0..M {
512 row_perm[i] = I::truncate(i);
513 }
514 for (i, t) in row_transpositions.iter().copied().enumerate() {
515 row_perm.as_mut().swap(i, t);
516 }
517 for i in 0..M {
518 row_perm_inv[row_perm[i].zx()] = I::truncate(i);
519 }
520 for j in 0..N {
521 col_perm[j] = I::truncate(j);
522 }
523 for (i, t) in col_transpositions.iter().copied().enumerate() {
524 col_perm.as_mut().swap(i, t);
525 }
526 for j in 0..N {
527 col_perm_inv[col_perm[j].zx()] = I::truncate(j);
528 }
529 unsafe {
530 (
531 FullPivLuInfo {
532 transposition_count: n_transpositions,
533 },
534 PermRef::new_unchecked(row_perm, row_perm_inv, M),
535 PermRef::new_unchecked(col_perm, col_perm_inv, N),
536 )
537 }
538}
539#[cfg(test)]
540mod tests {
541 use super::*;
542 use crate::stats::prelude::*;
543 use crate::utils::approx::*;
544 use crate::{Mat, assert, c64};
545 use dyn_stack::MemBuffer;
546 #[test]
547 fn test_flu() {
548 let rng = &mut StdRng::seed_from_u64(0);
549 for par in [Par::Seq, Par::rayon(8)] {
550 for m in [8, 16, 24, 32, 128, 255, 256, 257] {
551 let n = 8;
552 let approx_eq = CwiseMat(ApproxEq {
553 abs_tol: 1e-10,
554 rel_tol: 1e-10,
555 });
556 let A = CwiseMatDistribution {
557 nrows: m,
558 ncols: n,
559 dist: ComplexDistribution::new(
560 StandardNormal,
561 StandardNormal,
562 ),
563 }
564 .rand::<Mat<c64>>(rng);
565 let A = A.as_ref();
566 let mut LU = A.cloned();
567 let row_perm = &mut *vec![0usize; m];
568 let row_perm_inv = &mut *vec![0usize; m];
569 let col_perm = &mut *vec![0usize; n];
570 let col_perm_inv = &mut *vec![0usize; n];
571 let (_, p, q) = lu_in_place(
572 LU.as_mut(),
573 row_perm,
574 row_perm_inv,
575 col_perm,
576 col_perm_inv,
577 par,
578 MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
579 usize,
580 c64,
581 >(
582 n, n, par, default()
583 ))),
584 default(),
585 );
586 let mut L = LU.as_ref().cloned();
587 let mut U = LU.as_ref().cloned();
588 for j in 0..n {
589 for i in 0..j {
590 L[(i, j)] = c64::ZERO;
591 }
592 L[(j, j)] = c64::ONE;
593 }
594 for j in 0..n {
595 for i in j + 1..m {
596 U[(i, j)] = c64::ZERO;
597 }
598 }
599 let L = L.as_ref();
600 let U = U.as_ref();
601 let U = U.subrows(0, n);
602 assert!(p.inverse() * L * U * q ~ A);
603 }
604 for n in [16, 24, 32, 128, 255, 256, 257] {
605 let approx_eq = CwiseMat(ApproxEq {
606 abs_tol: 1e-10,
607 rel_tol: 1e-10,
608 });
609 let A = CwiseMatDistribution {
610 nrows: n,
611 ncols: n,
612 dist: StandardNormal,
613 }
614 .rand::<Mat<f64>>(rng);
615 let A = A.as_ref();
616 let mut LU = A.cloned();
617 let row_perm = &mut *vec![0usize; n];
618 let row_perm_inv = &mut *vec![0usize; n];
619 let col_perm = &mut *vec![0usize; n];
620 let col_perm_inv = &mut *vec![0usize; n];
621 let (_, p, q) = lu_in_place(
622 LU.as_mut(),
623 row_perm,
624 row_perm_inv,
625 col_perm,
626 col_perm_inv,
627 par,
628 MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
629 usize,
630 f64,
631 >(
632 n, n, par, default()
633 ))),
634 default(),
635 );
636 let mut L = LU.as_ref().cloned();
637 let mut U = LU.as_ref().cloned();
638 for j in 0..n {
639 for i in 0..j {
640 L[(i, j)] = 0.0;
641 }
642 L[(j, j)] = 1.0;
643 }
644 for j in 0..n {
645 for i in j + 1..n {
646 U[(i, j)] = 0.0;
647 }
648 }
649 let L = L.as_ref();
650 let U = U.as_ref();
651 assert!(p.inverse() * L * U * q ~ A);
652 }
653 }
654 }
655 #[test]
656 fn test_gh238() {
657 #[rustfmt::skip]
658 let A = [
659 [
660 -2.0,
661 0.0,
662 1.0,
663 0.0,
664 0.0,
665 0.0,
666 0.0,
667 0.0,
668 0.0,
669 0.0,
670 0.0,
671 0.0,
672 0.0,
673 0.0,
674 0.0,
675 0.0,
676 0.0,
677 0.0,
678 0.0,
679 0.0,
680 0.0,
681 0.0,
682 0.0,
683 0.0,
684 0.0,
685 0.0,
686 0.0,
687 0.0,
688 ],
689 [
690 0.0,
691 -2.0,
692 0.0,
693 1.0,
694 0.0,
695 0.0,
696 0.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 ],
719 [
720 0.0,
721 0.0,
722 0.0,
723 0.0,
724 -2.0,
725 0.0,
726 1.0,
727 0.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 ],
749 [
750 0.0,
751 0.0,
752 0.0,
753 0.0,
754 0.0,
755 -2.0,
756 0.0,
757 1.0,
758 0.0,
759 0.0,
760 0.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 ],
779 [
780 0.0,
781 0.0,
782 0.0,
783 0.0,
784 0.0,
785 0.0,
786 0.0,
787 0.0,
788 -2.0,
789 0.0,
790 1.0,
791 0.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 ],
809 [
810 0.0,
811 0.0,
812 0.0,
813 0.0,
814 0.0,
815 0.0,
816 0.0,
817 0.0,
818 0.0,
819 -2.0,
820 0.0,
821 1.0,
822 0.0,
823 0.0,
824 0.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 ],
839 [
840 0.0,
841 0.0,
842 0.0,
843 0.0,
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 -2.0,
853 0.0,
854 1.0,
855 0.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 ],
869 [
870 0.0,
871 0.0,
872 0.0,
873 0.0,
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 -2.0,
884 0.0,
885 1.0,
886 0.0,
887 0.0,
888 0.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 ],
899 [
900 0.0,
901 0.0,
902 0.0,
903 0.0,
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 -2.0,
917 0.0,
918 1.0,
919 0.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 ],
929 [
930 0.0,
931 0.0,
932 0.0,
933 0.0,
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 -2.0,
948 0.0,
949 1.0,
950 0.0,
951 0.0,
952 0.0,
953 0.0,
954 0.0,
955 0.0,
956 0.0,
957 0.0,
958 ],
959 [
960 0.0,
961 0.0,
962 0.0,
963 0.0,
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 -2.0,
981 0.0,
982 1.0,
983 0.0,
984 0.0,
985 0.0,
986 0.0,
987 0.0,
988 ],
989 [
990 0.0,
991 0.0,
992 0.0,
993 0.0,
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 -2.0,
1012 0.0,
1013 1.0,
1014 0.0,
1015 0.0,
1016 0.0,
1017 0.0,
1018 ],
1019 [
1020 0.0,
1021 0.0,
1022 0.0,
1023 0.0,
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 -2.0,
1045 0.0,
1046 1.0,
1047 0.0,
1048 ],
1049 [
1050 0.0,
1051 0.0,
1052 0.0,
1053 0.0,
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 -2.0,
1076 0.0,
1077 1.0,
1078 ],
1079 [
1080 0.0,
1081 0.0,
1082 0.0,
1083 0.0,
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 1.0,
1108 ],
1109 [
1110 0.0,
1111 0.0,
1112 0.0,
1113 1.0,
1114 0.0,
1115 0.0,
1116 1.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 ],
1139 [
1140 0.0,
1141 -1.0,
1142 0.0,
1143 0.0,
1144 1.0,
1145 0.0,
1146 0.0,
1147 0.0,
1148 0.0,
1149 0.0,
1150 0.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 ],
1169 [
1170 0.0,
1171 0.0,
1172 0.0,
1173 0.0,
1174 0.0,
1175 0.0,
1176 0.0,
1177 0.0,
1178 0.0,
1179 0.0,
1180 0.0,
1181 0.0,
1182 0.0,
1183 0.0,
1184 0.0,
1185 1.0,
1186 0.0,
1187 0.0,
1188 1.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 ],
1199 [
1200 0.0,
1201 0.0,
1202 0.0,
1203 0.0,
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 -1.0,
1214 0.0,
1215 0.0,
1216 1.0,
1217 0.0,
1218 0.0,
1219 0.0,
1220 0.0,
1221 0.0,
1222 0.0,
1223 0.0,
1224 0.0,
1225 0.0,
1226 0.0,
1227 0.0,
1228 ],
1229 [
1230 0.0,
1231 0.0,
1232 0.0,
1233 0.0,
1234 0.0,
1235 0.0,
1236 0.0,
1237 1.0,
1238 0.0,
1239 0.0,
1240 1.0,
1241 0.0,
1242 0.0,
1243 0.0,
1244 0.0,
1245 0.0,
1246 0.0,
1247 0.0,
1248 0.0,
1249 0.0,
1250 0.0,
1251 0.0,
1252 0.0,
1253 0.0,
1254 0.0,
1255 0.0,
1256 0.0,
1257 0.0,
1258 ],
1259 [
1260 0.0,
1261 0.0,
1262 0.0,
1263 0.0,
1264 0.0,
1265 1.0,
1266 0.0,
1267 0.0,
1268 -1.0,
1269 0.0,
1270 0.0,
1271 0.0,
1272 0.0,
1273 0.0,
1274 0.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 ],
1289 [
1290 0.0,
1291 0.0,
1292 0.0,
1293 0.0,
1294 0.0,
1295 0.0,
1296 0.0,
1297 0.0,
1298 0.0,
1299 0.0,
1300 0.0,
1301 1.0,
1302 0.0,
1303 0.0,
1304 1.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 ],
1319 [
1320 0.0,
1321 0.0,
1322 0.0,
1323 0.0,
1324 0.0,
1325 0.0,
1326 0.0,
1327 0.0,
1328 0.0,
1329 1.0,
1330 0.0,
1331 0.0,
1332 -1.0,
1333 0.0,
1334 0.0,
1335 0.0,
1336 0.0,
1337 0.0,
1338 0.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 ],
1349 [
1350 0.0,
1351 0.0,
1352 0.0,
1353 0.0,
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 0.0,
1364 0.0,
1365 0.0,
1366 0.0,
1367 0.0,
1368 0.0,
1369 0.0,
1370 0.0,
1371 0.0,
1372 0.0,
1373 1.0,
1374 0.0,
1375 0.0,
1376 1.0,
1377 0.0,
1378 ],
1379 [
1380 0.0,
1381 0.0,
1382 0.0,
1383 0.0,
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 -1.0,
1402 0.0,
1403 0.0,
1404 1.0,
1405 0.0,
1406 0.0,
1407 0.0,
1408 ],
1409 [
1410 0.0,
1411 0.0,
1412 1.0,
1413 0.0,
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 0.0,
1436 0.0,
1437 0.0,
1438 ],
1439 [
1440 0.0,
1441 0.0,
1442 0.0,
1443 0.0,
1444 0.0,
1445 0.0,
1446 0.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 1.0,
1460 0.0,
1461 0.0,
1462 1.0,
1463 0.0,
1464 0.0,
1465 0.0,
1466 0.0,
1467 0.0,
1468 ],
1469 [
1470 0.0,
1471 0.0,
1472 0.0,
1473 0.0,
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 1.0,
1488 0.0,
1489 0.0,
1490 -1.0,
1491 0.0,
1492 0.0,
1493 0.0,
1494 0.0,
1495 0.0,
1496 0.0,
1497 0.0,
1498 ],
1499 ];
1500 let A = MatRef::from_row_major_array(&A).to_owned();
1501 let b = [
1502 -0.0, -0.0, -2.0, -2.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0,
1503 -0.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1504 0.0, 0.0, 0.0, 0.0,
1505 ];
1506 let b = ColRef::from_slice(&b);
1507 let lu = A.full_piv_lu();
1508 assert!(lu.solve(b).is_all_finite());
1509 }
1510}