LU via Gauss Elimination with Complete Pivoting
(m: &mut ComplexMatrix)
| 2706 | |
| 2707 | /// LU via Gauss Elimination with Complete Pivoting |
| 2708 | fn gecp(m: &mut ComplexMatrix) -> (Vec<usize>, Vec<usize>) { |
| 2709 | let n = m.col; |
| 2710 | let mut r = vec![0usize; n - 1]; |
| 2711 | let mut s = vec![0usize; n - 1]; |
| 2712 | for k in 0..n - 1 { |
| 2713 | // Find pivot |
| 2714 | let (r_k, s_k) = match m.shape { |
| 2715 | Shape::Col => { |
| 2716 | let mut row_ics = 0usize; |
| 2717 | let mut col_ics = 0usize; |
| 2718 | let mut max_val = 0f64; |
| 2719 | for i in k..n { |
| 2720 | let c = m |
| 2721 | .col(i) |
| 2722 | .into_iter() |
| 2723 | .skip(k) |
| 2724 | .enumerate() |
| 2725 | .max_by(|x1, x2| x1.1.norm().partial_cmp(&x2.1.norm()).unwrap()) |
| 2726 | .unwrap(); |
| 2727 | let c_ics = c.0 + k; |
| 2728 | let c_val = c.1.norm(); |
| 2729 | if c_val > max_val { |
| 2730 | row_ics = c_ics; |
| 2731 | col_ics = i; |
| 2732 | max_val = c_val; |
| 2733 | } |
| 2734 | } |
| 2735 | (row_ics, col_ics) |
| 2736 | } |
| 2737 | Shape::Row => { |
| 2738 | let mut row_ics = 0usize; |
| 2739 | let mut col_ics = 0usize; |
| 2740 | let mut max_val = 0f64; |
| 2741 | for i in k..n { |
| 2742 | let c = m |
| 2743 | .row(i) |
| 2744 | .into_iter() |
| 2745 | .skip(k) |
| 2746 | .enumerate() |
| 2747 | .max_by(|x1, x2| x1.1.norm().partial_cmp(&x2.1.norm()).unwrap()) |
| 2748 | .unwrap(); |
| 2749 | let c_ics = c.0 + k; |
| 2750 | let c_val = c.1.norm(); |
| 2751 | if c_val > max_val { |
| 2752 | col_ics = c_ics; |
| 2753 | row_ics = i; |
| 2754 | max_val = c_val; |
| 2755 | } |
| 2756 | } |
| 2757 | (row_ics, col_ics) |
| 2758 | } |
| 2759 | }; |
| 2760 | r[k] = r_k; |
| 2761 | s[k] = s_k; |
| 2762 | |
| 2763 | // Interchange rows |
| 2764 | for j in k..n { |
| 2765 | unsafe { |