fact: Fixpoint cormen_lup {m n} := match m, n return 'M_(m.+1,n.+1) -> 'S_m.+1 * 'M_(m.+1,m.+1) * 'M_(m.+1,n.+1) with | p.+1, _.+1 => fun (A : 'M_(1 + (1 + p), 1 + _)) => let k := odflt 0 (find_pivot A) in let A1 : 'M_(1 + _, 1 + _) := xrow 0 k A in let P1 : 'S_(1 + (1 + p)) := tperm 0 k in let Schur := ((fun_of_matrix A k 0)^-1 *: dlsubmx A1) *m ursubmx A1 in let: (P2, L2, U2) := cormen_lup (drsubmx A1 - Schur) in let P := (lift0_perm P2 * P1)%g in let pA1 := row_perm P2 (dlsubmx A1) in let L := block_mx 1%:M (const_mx 0) ((fun_of_matrix A k 0)^-1 *: pA1) L2 in let U := block_mx (ulsubmx A1) (ursubmx A1) (const_mx 0) U2 in (P, L, U) | _, _ => fun A => (1%g, 1%:M, A) end.
imports: mathcomp(ssreflect ssrbool ssrfun eqtype ssrnat seq choice fintype) mathcomp(div finfun bigop prime binomial ssralg finset fingroup finalg) mathcomp(mxalgebra perm zmodp matrix ssrint)