fn gecp(m: &mut Matrix) -> (Vec<usize>, Vec<usize>) {
    let n = m.col;
    let mut r = vec![0usize; n - 1];
    let mut s = vec![0usize; n - 1];
    for k in 0..n - 1 {
        // Find pivot
        let (r_k, s_k) = match m.shape {
            Col => {
                let mut row_ics = 0usize;
                let mut col_ics = 0usize;
                let mut max_val = 0f64;
                for i in k..n {
                    let c = m
                        .col(i)
                        .into_iter()
                        .skip(k)
                        .enumerate()
                        .max_by(|x1, x2| x1.1.abs().partial_cmp(&x2.1.abs()).unwrap())
                        .unwrap();
                    let c_ics = c.0 + k;
                    let c_val = c.1.abs();
                    if c_val > max_val {
                        row_ics = c_ics;
                        col_ics = i;
                        max_val = c_val;
                    }
                }
                (row_ics, col_ics)
            }
            Row => {
                let mut row_ics = 0usize;
                let mut col_ics = 0usize;
                let mut max_val = 0f64;
                for i in k..n {
                    let c = m
                        .row(i)
                        .into_iter()
                        .skip(k)
                        .enumerate()
                        .max_by(|x1, x2| x1.1.abs().partial_cmp(&x2.1.abs()).unwrap())
                        .unwrap();
                    let c_ics = c.0 + k;
                    let c_val = c.1.abs();
                    if c_val > max_val {
                        col_ics = c_ics;
                        row_ics = i;
                        max_val = c_val;
                    }
                }
                (row_ics, col_ics)
            }
        };
        r[k] = r_k;
        s[k] = s_k;

        // Interchange rows
        for j in k..n {
            unsafe {
                std::ptr::swap(&mut m[(k, j)], &mut m[(r_k, j)]);
            }
        }

        // Interchange cols
        for i in 0..n {
            unsafe {
                std::ptr::swap(&mut m[(i, k)], &mut m[(i, s_k)]);
            }
        }

        // Form the multipliers
        for i in k + 1..n {
            m[(i, k)] = -m[(i, k)] / m[(k, k)];
            for j in k + 1..n {
                m[(i, j)] += m[(i, k)] * m[(k, j)];
            }
        }
    }
    (r, s)
}
LU via Gauss Elimination with Complete Pivoting
src/structure/matrix.rs