fn lu(&self) -> PQLU<ComplexMatrix> {
        assert_eq!(self.col, self.row);
        let n = self.row;
        let len: usize = n * n;

        let mut l = ceye(n);
        let mut u = cmatrix(vec![Complex::zero(); len], n, n, self.shape);

        let mut temp = self.clone();
        let (p, q) = gecp(&mut temp);
        for i in 0..n {
            for j in 0..i {
                // Inverse multiplier
                l[(i, j)] = -temp[(i, j)];
            }
            for j in i..n {
                u[(i, j)] = temp[(i, j)];
            }
        }
        // Pivoting L
        for i in 0..n - 1 {
            unsafe {
                let l_i = l.col_mut(i);
                for j in i + 1..l.col - 1 {
                    let dst = p[j];
                    std::ptr::swap(l_i[j], l_i[dst]);
                }
            }
        }
        PQLU { p, q, l, u }
    }
LU Decomposition Implements (Complete Pivot)

# Description
It use complete pivoting LU decomposition.
You can get two permutations, and LU matrices.

# Caution
It returns `Option<PQLU>` - You should unwrap to obtain real value.
`PQLU` has four field - `p`, `q`, `l`, `u`.
`p`, `q` are permutations.
`l`, `u` are matrices.

# Examples
```
#[macro_use]
use peroxide::fuga::*;
use peroxide::complex::matrix::*;

fn main() {
let a = cmatrix(vec![
C64::new(1f64, 1f64),
C64::new(2f64, 2f64),
C64::new(3f64, 3f64),
C64::new(4f64, 4f64)
],
2, 2, Row
);

let l_exp = cmatrix(vec![
C64::new(1f64, 0f64),
C64::new(0f64, 0f64),
C64::new(0.5f64, -0.0f64),
C64::new(1f64, 0f64)
],
2, 2, Row
);

let u_exp = cmatrix(vec![
C64::new(4f64, 4f64),
C64::new(3f64, 3f64),
C64::new(0f64, 0f64),
C64::new(-0.5f64, -0.5f64)
],
2, 2, Row
);
let pqlu = a.lu();
let (p,q,l,u) = (pqlu.p, pqlu.q, pqlu.l, pqlu.u);
assert_eq!(p, vec![1]); // swap 0 & 1 (Row)
assert_eq!(q, vec![1]); // swap 0 & 1 (Col)
assert_eq!(l, l_exp);
assert_eq!(u, u_exp);
}
```
src/complex/matrix.rs