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();
488 if mat.row_stride().unsigned_abs() != 1
489 && M >= N
490 && N > 1
491 && linalg::row_major_copy_pays(M, N)
492 {
493 let mut work = alloc::vec![zero::<T>(); M * N];
494 let mut work = MatMut::from_column_major_slice_mut(&mut work, M, N);
495 linalg::copy_across_layouts(work.rb_mut(), mat.rb(), par);
496 let ret = lu_in_place(
497 work.rb_mut(),
498 row_perm,
499 row_perm_inv,
500 col_perm,
501 col_perm_inv,
502 par,
503 stack,
504 params,
505 );
506 linalg::copy_across_layouts(mat, work.rb(), par);
507 return ret;
508 }
509 let size = Ord::min(M, N);
510 let (mut row_transpositions, stack) = stack.make_with(size, |_| 0);
511 let row_transpositions = row_transpositions.as_mut();
512 let (mut col_transpositions, _) = stack.make_with(size, |_| 0);
513 let col_transpositions = col_transpositions.as_mut();
514 let n_transpositions = if mat.row_stride().abs() < mat.col_stride().abs() {
515 lu_in_place_unblocked(
516 mat,
517 row_transpositions,
518 col_transpositions,
519 par,
520 false,
521 params,
522 )
523 } else {
524 lu_in_place_unblocked(
525 mat.transpose_mut(),
526 col_transpositions,
527 row_transpositions,
528 par,
529 true,
530 params,
531 )
532 };
533 for i in 0..M {
534 row_perm[i] = I::truncate(i);
535 }
536 for (i, t) in row_transpositions.iter().copied().enumerate() {
537 row_perm.as_mut().swap(i, t);
538 }
539 for i in 0..M {
540 row_perm_inv[row_perm[i].zx()] = I::truncate(i);
541 }
542 for j in 0..N {
543 col_perm[j] = I::truncate(j);
544 }
545 for (i, t) in col_transpositions.iter().copied().enumerate() {
546 col_perm.as_mut().swap(i, t);
547 }
548 for j in 0..N {
549 col_perm_inv[col_perm[j].zx()] = I::truncate(j);
550 }
551 unsafe {
552 (
553 FullPivLuInfo {
554 transposition_count: n_transpositions,
555 },
556 PermRef::new_unchecked(row_perm, row_perm_inv, M),
557 PermRef::new_unchecked(col_perm, col_perm_inv, N),
558 )
559 }
560}
561#[cfg(test)]
562mod tests {
563 use super::*;
564 use crate::stats::prelude::*;
565 use crate::utils::approx::*;
566 use crate::{Mat, assert, c64};
567 use dyn_stack::MemBuffer;
568 #[test]
569 fn test_flu() {
570 let rng = &mut StdRng::seed_from_u64(0);
571 for par in [Par::Seq, Par::rayon(8)] {
572 for m in [8, 16, 24, 32, 128, 255, 256, 257] {
573 let n = 8;
574 let approx_eq = CwiseMat(ApproxEq {
575 abs_tol: 1e-10,
576 rel_tol: 1e-10,
577 });
578 let A = CwiseMatDistribution {
579 nrows: m,
580 ncols: n,
581 dist: ComplexDistribution::new(
582 StandardNormal,
583 StandardNormal,
584 ),
585 }
586 .rand::<Mat<c64>>(rng);
587 let A = A.as_ref();
588 let mut LU = A.cloned();
589 let row_perm = &mut *vec![0usize; m];
590 let row_perm_inv = &mut *vec![0usize; m];
591 let col_perm = &mut *vec![0usize; n];
592 let col_perm_inv = &mut *vec![0usize; n];
593 let (_, p, q) = lu_in_place(
594 LU.as_mut(),
595 row_perm,
596 row_perm_inv,
597 col_perm,
598 col_perm_inv,
599 par,
600 MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
601 usize,
602 c64,
603 >(
604 n, n, par, default()
605 ))),
606 default(),
607 );
608 let mut L = LU.as_ref().cloned();
609 let mut U = LU.as_ref().cloned();
610 for j in 0..n {
611 for i in 0..j {
612 L[(i, j)] = c64::ZERO;
613 }
614 L[(j, j)] = c64::ONE;
615 }
616 for j in 0..n {
617 for i in j + 1..m {
618 U[(i, j)] = c64::ZERO;
619 }
620 }
621 let L = L.as_ref();
622 let U = U.as_ref();
623 let U = U.subrows(0, n);
624 assert!(p.inverse() * L * U * q ~ A);
625 }
626 for n in [16, 24, 32, 128, 255, 256, 257] {
627 let approx_eq = CwiseMat(ApproxEq {
628 abs_tol: 1e-10,
629 rel_tol: 1e-10,
630 });
631 let A = CwiseMatDistribution {
632 nrows: n,
633 ncols: n,
634 dist: StandardNormal,
635 }
636 .rand::<Mat<f64>>(rng);
637 let A = A.as_ref();
638 let mut LU = A.cloned();
639 let row_perm = &mut *vec![0usize; n];
640 let row_perm_inv = &mut *vec![0usize; n];
641 let col_perm = &mut *vec![0usize; n];
642 let col_perm_inv = &mut *vec![0usize; n];
643 let (_, p, q) = lu_in_place(
644 LU.as_mut(),
645 row_perm,
646 row_perm_inv,
647 col_perm,
648 col_perm_inv,
649 par,
650 MemStack::new(&mut MemBuffer::new(lu_in_place_scratch::<
651 usize,
652 f64,
653 >(
654 n, n, par, default()
655 ))),
656 default(),
657 );
658 let mut L = LU.as_ref().cloned();
659 let mut U = LU.as_ref().cloned();
660 for j in 0..n {
661 for i in 0..j {
662 L[(i, j)] = 0.0;
663 }
664 L[(j, j)] = 1.0;
665 }
666 for j in 0..n {
667 for i in j + 1..n {
668 U[(i, j)] = 0.0;
669 }
670 }
671 let L = L.as_ref();
672 let U = U.as_ref();
673 assert!(p.inverse() * L * U * q ~ A);
674 }
675 }
676 }
677 #[test]
678 fn test_gh238() {
679 #[rustfmt::skip]
680 let A = [
681 [
682 -2.0,
683 0.0,
684 1.0,
685 0.0,
686 0.0,
687 0.0,
688 0.0,
689 0.0,
690 0.0,
691 0.0,
692 0.0,
693 0.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 ],
711 [
712 0.0,
713 -2.0,
714 0.0,
715 1.0,
716 0.0,
717 0.0,
718 0.0,
719 0.0,
720 0.0,
721 0.0,
722 0.0,
723 0.0,
724 0.0,
725 0.0,
726 0.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 ],
741 [
742 0.0,
743 0.0,
744 0.0,
745 0.0,
746 -2.0,
747 0.0,
748 1.0,
749 0.0,
750 0.0,
751 0.0,
752 0.0,
753 0.0,
754 0.0,
755 0.0,
756 0.0,
757 0.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 ],
771 [
772 0.0,
773 0.0,
774 0.0,
775 0.0,
776 0.0,
777 -2.0,
778 0.0,
779 1.0,
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 0.0,
789 0.0,
790 0.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 ],
801 [
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 -2.0,
811 0.0,
812 1.0,
813 0.0,
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 0.0,
823 0.0,
824 0.0,
825 0.0,
826 0.0,
827 0.0,
828 0.0,
829 0.0,
830 ],
831 [
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 -2.0,
842 0.0,
843 1.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 0.0,
853 0.0,
854 0.0,
855 0.0,
856 0.0,
857 0.0,
858 0.0,
859 0.0,
860 ],
861 [
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 0.0,
873 0.0,
874 -2.0,
875 0.0,
876 1.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 0.0,
887 0.0,
888 0.0,
889 0.0,
890 ],
891 [
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 0.0,
903 0.0,
904 0.0,
905 -2.0,
906 0.0,
907 1.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 0.0,
918 0.0,
919 0.0,
920 ],
921 [
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 0.0,
933 0.0,
934 0.0,
935 0.0,
936 0.0,
937 0.0,
938 -2.0,
939 0.0,
940 1.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 ],
951 [
952 0.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 0.0,
963 0.0,
964 0.0,
965 0.0,
966 0.0,
967 0.0,
968 0.0,
969 -2.0,
970 0.0,
971 1.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 ],
981 [
982 0.0,
983 0.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 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 -2.0,
1003 0.0,
1004 1.0,
1005 0.0,
1006 0.0,
1007 0.0,
1008 0.0,
1009 0.0,
1010 ],
1011 [
1012 0.0,
1013 0.0,
1014 0.0,
1015 0.0,
1016 0.0,
1017 0.0,
1018 0.0,
1019 0.0,
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 -2.0,
1034 0.0,
1035 1.0,
1036 0.0,
1037 0.0,
1038 0.0,
1039 0.0,
1040 ],
1041 [
1042 0.0,
1043 0.0,
1044 0.0,
1045 0.0,
1046 0.0,
1047 0.0,
1048 0.0,
1049 0.0,
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 -2.0,
1067 0.0,
1068 1.0,
1069 0.0,
1070 ],
1071 [
1072 0.0,
1073 0.0,
1074 0.0,
1075 0.0,
1076 0.0,
1077 0.0,
1078 0.0,
1079 0.0,
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 -2.0,
1098 0.0,
1099 1.0,
1100 ],
1101 [
1102 0.0,
1103 0.0,
1104 0.0,
1105 0.0,
1106 0.0,
1107 0.0,
1108 0.0,
1109 0.0,
1110 0.0,
1111 0.0,
1112 0.0,
1113 0.0,
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 1.0,
1130 ],
1131 [
1132 0.0,
1133 0.0,
1134 0.0,
1135 1.0,
1136 0.0,
1137 0.0,
1138 1.0,
1139 0.0,
1140 0.0,
1141 0.0,
1142 0.0,
1143 0.0,
1144 0.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 ],
1161 [
1162 0.0,
1163 -1.0,
1164 0.0,
1165 0.0,
1166 1.0,
1167 0.0,
1168 0.0,
1169 0.0,
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 0.0,
1186 0.0,
1187 0.0,
1188 0.0,
1189 0.0,
1190 ],
1191 [
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 0.0,
1203 0.0,
1204 0.0,
1205 0.0,
1206 0.0,
1207 1.0,
1208 0.0,
1209 0.0,
1210 1.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 0.0,
1220 ],
1221 [
1222 0.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 0.0,
1233 0.0,
1234 0.0,
1235 -1.0,
1236 0.0,
1237 0.0,
1238 1.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 0.0,
1248 0.0,
1249 0.0,
1250 ],
1251 [
1252 0.0,
1253 0.0,
1254 0.0,
1255 0.0,
1256 0.0,
1257 0.0,
1258 0.0,
1259 1.0,
1260 0.0,
1261 0.0,
1262 1.0,
1263 0.0,
1264 0.0,
1265 0.0,
1266 0.0,
1267 0.0,
1268 0.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 ],
1281 [
1282 0.0,
1283 0.0,
1284 0.0,
1285 0.0,
1286 0.0,
1287 1.0,
1288 0.0,
1289 0.0,
1290 -1.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 0.0,
1302 0.0,
1303 0.0,
1304 0.0,
1305 0.0,
1306 0.0,
1307 0.0,
1308 0.0,
1309 0.0,
1310 ],
1311 [
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 0.0,
1323 1.0,
1324 0.0,
1325 0.0,
1326 1.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 0.0,
1336 0.0,
1337 0.0,
1338 0.0,
1339 0.0,
1340 ],
1341 [
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 1.0,
1352 0.0,
1353 0.0,
1354 -1.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 ],
1371 [
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 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 1.0,
1396 0.0,
1397 0.0,
1398 1.0,
1399 0.0,
1400 ],
1401 [
1402 0.0,
1403 0.0,
1404 0.0,
1405 0.0,
1406 0.0,
1407 0.0,
1408 0.0,
1409 0.0,
1410 0.0,
1411 0.0,
1412 0.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 -1.0,
1424 0.0,
1425 0.0,
1426 1.0,
1427 0.0,
1428 0.0,
1429 0.0,
1430 ],
1431 [
1432 0.0,
1433 0.0,
1434 1.0,
1435 0.0,
1436 0.0,
1437 0.0,
1438 0.0,
1439 0.0,
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 0.0,
1460 ],
1461 [
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 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 1.0,
1482 0.0,
1483 0.0,
1484 1.0,
1485 0.0,
1486 0.0,
1487 0.0,
1488 0.0,
1489 0.0,
1490 ],
1491 [
1492 0.0,
1493 0.0,
1494 0.0,
1495 0.0,
1496 0.0,
1497 0.0,
1498 0.0,
1499 0.0,
1500 0.0,
1501 0.0,
1502 0.0,
1503 0.0,
1504 0.0,
1505 0.0,
1506 0.0,
1507 0.0,
1508 0.0,
1509 1.0,
1510 0.0,
1511 0.0,
1512 -1.0,
1513 0.0,
1514 0.0,
1515 0.0,
1516 0.0,
1517 0.0,
1518 0.0,
1519 0.0,
1520 ],
1521 ];
1522 let A = MatRef::from_row_major_array(&A).to_owned();
1523 let b = [
1524 -0.0, -0.0, -2.0, -2.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0, -0.0,
1525 -0.0, -0.0, -0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
1526 0.0, 0.0, 0.0, 0.0,
1527 ];
1528 let b = ColRef::from_slice(&b);
1529 let lu = A.full_piv_lu();
1530 assert!(lu.solve(b).is_all_finite());
1531 }
1532}