Library mathcomp.algebra.matrix

(* (c) Copyright 2006-2016 Microsoft Corporation and Inria.                  
 Distributed under the terms of CeCILL-B.                                  *)

From HB Require Import structures.
From mathcomp Require Import ssreflect ssrbool ssrfun eqtype ssrnat seq choice.
From mathcomp Require Import fintype finfun bigop finset fingroup perm order.
From mathcomp Require Import div prime binomial ssralg countalg finalg zmodp.

Basic concrete linear algebra : definition of type for matrices, and all basic matrix operations including determinant, trace and support for block decomposition. Matrices are represented by a row-major list of their coefficients but this implementation is hidden by three levels of wrappers (Matrix/Finfun/Tuple) so the matrix type should be treated as abstract and handled using only the operations described below: 'M[R](m, n) == the type of m rows by n columns matrices with 'M(m, n) coefficients in R; the [R] is optional and is usually omitted. 'M[R]_n, 'M_n == the type of n x n square matrices. 'rV[R]_n, 'rV_n == the type of 1 x n row vectors. 'cV[R]_n, 'cV_n == the type of n x 1 column vectors. \matrix(i < m, j < n) Expr(i, j) == the m x n matrix with general coefficient Expr(i, j), with i : 'I_m and j : 'I_n. the < m bound can be omitted if it is equal to n, though usually both bounds are omitted as they can be inferred from the context. \row(j < n) Expr(j), \col(i < m) Expr(i) the row / column vectors with general term Expr; the parentheses can be omitted along with the bound. \matrix(i < m) RowExpr(i) == the m x n matrix with row i given by RowExpr(i) : 'rV_n. A i j == the coefficient of matrix A : 'M(m, n) in column j of row i, where i : 'I_m, and j : 'I_n (via the coercion fun_of_matrix : matrix >-> Funclass). const_mx a == the constant matrix whose entries are all a (dimensions should be determined by context). map_mx f A == the pointwise image of A by f, i.e., the matrix Af congruent to A with Af i j = f (A i j) for all i and j. A^T == the matrix transpose of A. row i A == the i'th row of A (this is a row vector). col j A == the j'th column of A (a column vector). row' i A == A with the i'th row spliced out. col' i A == A with the j'th column spliced out. xrow i1 i2 A == A with rows i1 and i2 interchanged. xcol j1 j2 A == A with columns j1 and j2 interchanged. row_perm s A == A : 'M(m, n) with rows permuted by s : 'S_m. col_perm s A == A : 'M(m, n) with columns permuted by s : 'S_n. row_mx Al Ar == the row block matrix <Al Ar> obtained by contatenating two matrices Al and Ar of the same height. col_mx Au Ad == the column block matrix / Au \ (Au and Ad must have the same width). \ Ad / block_mx Aul Aur Adl Adr == the block matrix / Aul Aur \ \ Adl Adr / \mxblock(i < m, j < n) B i j == the block matrix of type 'M(\sum_i p i, \sum_j q j) / (B 0 0) ⋯ (B 0 j) ⋯ (B 0 n) \ | ... ... ... | | (B i 0) ⋯ (B i j) ⋯ (B i n) | | ... ... ... | \ (B m 0) ⋯ (B m j) ⋯ (B m n) / where each block (B i j) has type 'M(p i, q j). \mxdiag(i < n) B i == the block square matrix of type 'M(\sum_i p i) / (B 0) 0 \ | ... ... | | 0 (B i) 0 | | ... ... | \ 0 (B n) / where each block (B i) has type 'M(p i). \mxrow(j < n) B j == the block matrix of type 'M(m, \sum_j q j). < (B 0) ... (B n) > where each block (B j) has type 'M(m, q j). \mxcol(i < m) B i == the block matrix of type 'M(\sum_i p i, n) / (B 0) \ | ... | \ (B m) / where each block (B i) has type 'M(p i, n). [l|r]submx A == the left/right submatrices of a row block matrix A. Note that the type of A, 'M(m, n1 + n2) indicates how A should be decomposed. [u|d]submx A == the up/down submatrices of a column block matrix A. [u|d] [l|r]submx A == the upper left, etc submatrices of a block matrix A. submxblock A i j == the block submatrix of type 'M(p i, q j) of A. The type of A, 'M(\sum_i p i, \sum_i q i) indicates how A should be decomposed. There is no analogous for mxdiag since one can use submxblock A i i to extract a diagonal block. submxrow A j == the submatrix of type 'M(m, q j) of A. The type of A, 'M(m, \sum_j q j) indicates how A should be decomposed. submxrow A j == the submatrix of type 'M(p i, n) of A. The type of A, 'M(\sum_i p i, n) indicates how A should be decomposed. mxsub f g A == generic reordered submatrix, given by functions f and g which specify which subset of rows and columns to take and how to reorder them, e.g. picking f and g to be increasing yields traditional submatrices. := \matrix(i, j) A (f i) (g i) rowsub f A := mxsub f id A colsub g A := mxsub id g A castmx eq_mn A == A : 'M(m, n) cast to 'M(m', n') using the equation pair eq_mn : (m = m') * (n = n'). This is the usual workaround for the syntactic limitations of dependent types in Coq, and can be used to introduce a block decomposition. It simplifies to A when eq_mn is the pair (erefl m, erefl n) (using rewrite /castmx /=). conform_mx B A == A if A and B have the same dimensions, else B. mxvec A == a row vector of width m * n holding all the entries of the m x n matrix A. mxvec_index i j == the index of A i j in mxvec A. vec_mx v == the inverse of mxvec, reshaping a vector of width m * n back into into an m x n rectangular matrix. In 'M[R](m, n), R can be any type, but 'M[R](m, n) inherits the eqType, choiceType, countType, finType, zmodType structures of R; 'M[R](m, n) also has a natural lmodType R structure when R has a ringType structure. Because the type of matrices specifies their dimension, only non-trivial square matrices (of type 'M[R]_n.+1) can inherit the ring structure of R; indeed they then have an algebra structure (lalgType R, or algType R if R is a comRingType, or even unitAlgType if R is a comUnitRingType). We thus provide separate syntax for the general matrix multiplication, and other operations for matrices over a ringType R: A *m B == the matrix product of A and B; the width of A must be equal to the height of B. a%:M == the scalar matrix with a's on the main diagonal; in particular 1%:M denotes the identity matrix, and is is equal to 1%R when n is of the form n'.+1 (e.g., n >= 1). is_scalar_mx A <=> A is a scalar matrix (A = a%:M for some A). diag_mx d == the diagonal matrix whose main diagonal is d : 'rV_n. is_diag_mx A <=> A is a diagonal matrix: forall i j, i != j -> A i j = 0 is_trig_mx A <=> A is a triangular matrix: forall i j, i < j -> A i j = 0 delta_mx i j == the matrix with a 1 in row i, column j and 0 elsewhere. pid_mx r == the partial identity matrix with 1s only on the r first coefficients of the main diagonal; the dimensions of pid_mx r are determined by the context, and pid_mx r can be rectangular. copid_mx r == the complement to 1%:M of pid_mx r: a square diagonal matrix with 1s on all but the first r coefficients on its main diagonal. perm_mx s == the n x n permutation matrix for s : 'S_n. tperm_mx i1 i2 == the permutation matrix that exchanges i1 i2 : 'I_n. is_perm_mx A == A is a permutation matrix. lift0_mx A == the 1 + n square matrix block_mx 1 0 0 A when A : 'M_n. \tr A == the trace of a square matrix A. \det A == the determinant of A, using the Leibnitz formula. cofactor i j A == the i, j cofactor of A (the signed i, j minor of A), \adj A == the adjugate matrix of A (\adj A i j = cofactor j i A). A \in unitmx == A is invertible (R must be a comUnitRingType). invmx A == the inverse matrix of A if A \in unitmx A, otherwise A. A \is a mxOver S == the matrix A has its coefficients in S. comm_mx A B := A *m B = B *m A comm_mxb A B := A *m B == B *m A all_comm_mx As fs := all2rel comm_mxb fs The following operations provide a correspondence between linear functions and matrices: lin1_mx f == the m x n matrix that emulates via right product a (linear) function f : 'rV_m -> 'rV_n on ROW VECTORS lin_mx f == the (m1 * n1) x (m2 * n2) matrix that emulates, via the right multiplication on the mxvec encodings, a linear function f : 'M(m1, n1) -> 'M(m2, n2) lin_mul_row u := lin1_mx (mulmx u \o vec_mx) (applies a row-encoded function to the row-vector u). mulmx A == partially applied matrix multiplication (mulmx A B is displayed as A *m B), with, for A : 'M(m, n), a canonical {linear 'M(n, p) -> 'M(m, p}} structure. mulmxr A == self-simplifying right-hand matrix multiplication, i.e., mulmxr A B simplifies to B *m A, with, for A : 'M(n, p), a canonical {linear 'M(m, n) -> 'M(m, p}} structure. lin_mulmx A := lin_mx (mulmx A). lin_mulmxr A := lin_mx (mulmxr A). We also extend any finType structure of R to 'M[R](m, n), and define: {'GL_n[R]} == the finGroupType of units of 'M[R]_n.-1.+1. 'GL_n[R] == the general linear group of all matrices in {'GL_n(R)}. 'GL_n(p) == 'GL_n['F_p], the general linear group of a prime field. GLval u == the coercion of u : {'GL_n(R)} to a matrix. In addition to the lemmas relevant to these definitions, this file also proves several classic results, including :
  • The determinant is a multilinear alternate form.
  • The Laplace determinant expansion formulas: expand_det [row|col].
  • The Cramer rule : mul_mx_adj & mul_adj_mx.
Finally, as an example of the use of block products, we program and prove the correctness of a classical linear algebra algorithm: cormen_lup A == the triangular decomposition (L, U, P) of a nontrivial square matrix A into a lower triagular matrix L with 1s on the main diagonal, an upper matrix U, and a permutation matrix P, such that P * A = L * U. This is example only; we use a different, more precise algorithm to develop the theory of matrix ranks and row spaces in mxalgebra.v

Set Implicit Arguments.

Import GroupScope.
Import GRing.Theory.
Local Open Scope ring_scope.

Reserved Notation "''M_' n" (at level 8, n at level 2, format "''M_' n").
Reserved Notation "''rV_' n" (at level 8, n at level 2, format "''rV_' n").
Reserved Notation "''cV_' n" (at level 8, n at level 2, format "''cV_' n").
Reserved Notation "''M_' ( n )" (at level 8). (* only parsing *)
Reserved Notation "''M_' ( m , n )" (at level 8, format "''M_' ( m , n )").
Reserved Notation "''M[' R ]_ n" (at level 8, n at level 2). (* only parsing *)
Reserved Notation "''rV[' R ]_ n" (at level 8, n at level 2). (* only parsing *)
Reserved Notation "''cV[' R ]_ n" (at level 8, n at level 2). (* only parsing *)
Reserved Notation "''M[' R ]_ ( n )" (at level 8). (* only parsing *)
Reserved Notation "''M[' R ]_ ( m , n )" (at level 8). (* only parsing *)

Reserved Notation "\matrix_ i E"
  (at level 36, E at level 36, i at level 2,
   format "\matrix_ i E").
Reserved Notation "\matrix_ ( i < n ) E"
  (at level 36, E at level 36, i, n at level 50). (* only parsing *)
Reserved Notation "\matrix_ ( i , j ) E"
  (at level 36, E at level 36, i, j at level 50,
   format "\matrix_ ( i , j ) E").
Reserved Notation "\matrix[ k ]_ ( i , j ) E"
  (at level 36, E at level 36, i, j at level 50,
   format "\matrix[ k ]_ ( i , j ) E").
Reserved Notation "\matrix_ ( i < m , j < n ) E"
  (at level 36, E at level 36, i, m, j, n at level 50). (* only parsing *)
Reserved Notation "\matrix_ ( i , j < n ) E"
  (at level 36, E at level 36, i, j, n at level 50). (* only parsing *)
Reserved Notation "\row_ j E"
  (at level 36, E at level 36, j at level 2,
   format "\row_ j E").
Reserved Notation "\row_ ( j < n ) E"
  (at level 36, E at level 36, j, n at level 50). (* only parsing *)
Reserved Notation "\col_ j E"
  (at level 36, E at level 36, j at level 2,
   format "\col_ j E").
Reserved Notation "\col_ ( j < n ) E"
  (at level 36, E at level 36, j, n at level 50). (* only parsing *)
Reserved Notation "\mxblock_ ( i , j ) E"
  (at level 36, E at level 36, i, j at level 50,
   format "\mxblock_ ( i , j ) E").
Reserved Notation "\mxblock_ ( i < m , j < n ) E"
  (at level 36, E at level 36, i, m, j, n at level 50). (* only parsing *)
Reserved Notation "\mxblock_ ( i , j < n ) E"
  (at level 36, E at level 36, i, j, n at level 50). (* only parsing *)
Reserved Notation "\mxrow_ j E"
  (at level 36, E at level 36, j at level 2,
   format "\mxrow_ j E").
Reserved Notation "\mxrow_ ( j < n ) E"
  (at level 36, E at level 36, j, n at level 50). (* only parsing *)
Reserved Notation "\mxcol_ j E"
  (at level 36, E at level 36, j at level 2,
   format "\mxcol_ j E").
Reserved Notation "\mxcol_ ( j < n ) E"
  (at level 36, E at level 36, j, n at level 50). (* only parsing *)
Reserved Notation "\mxdiag_ j E"
  (at level 36, E at level 36, j at level 2,
   format "\mxdiag_ j E").
Reserved Notation "\mxdiag_ ( j < n ) E"
  (at level 36, E at level 36, j, n at level 50). (* only parsing *)

Reserved Notation "x %:M" (at level 8, format "x %:M").
Reserved Notation "A *m B" (at level 40, left associativity, format "A *m B").
Reserved Notation "A ^T" (at level 8, format "A ^T").
Reserved Notation "\tr A" (at level 10, A at level 8, format "\tr A").
Reserved Notation "\det A" (at level 10, A at level 8, format "\det A").
Reserved Notation "\adj A" (at level 10, A at level 8, format "\adj A").


Type Definition*********************************

Section MatrixDef.

Variable R : Type.
Variables m n : nat.

Basic linear algebra (matrices). We use dependent types (ordinals) for the indices so that ranges are mostly inferred automatically

Variant matrix : predArgType := Matrix of {ffun 'I_m × 'I_n R}.

Definition mx_val A := let: Matrix g := A in g.


Definition fun_of_matrix A (i : 'I_m) (j : 'I_n) := mx_val A (i, j).

Coercion fun_of_matrix : matrix >-> Funclass.

End MatrixDef.

Fact matrix_key : unit.

Canonical matrix_unlockable := Unlockable matrix_of_fun.unlock.

Section MatrixDef2.

Variable R : Type.
Variables m n : nat.
Implicit Type F : 'I_m 'I_n R.

Lemma mxE k F : matrix_of_fun k F =2 F.

Lemma matrixP (A B : matrix R m n) : A =2 B A = B.

Lemma eq_mx k F1 F2 : (F1 =2 F2) matrix_of_fun k F1 = matrix_of_fun k F2.

End MatrixDef2.

Arguments eq_mx {R m n k} [F1] F2 eq_F12.

Bind Scope ring_scope with matrix.

Notation "''M[' R ]_ ( m , n )" := (matrix R m n) (only parsing): type_scope.
Notation "''rV[' R ]_ n" := 'M[R]_(1, n) (only parsing) : type_scope.
Notation "''cV[' R ]_ n" := 'M[R]_(n, 1) (only parsing) : type_scope.
Notation "''M[' R ]_ n" := 'M[R]_(n, n) (only parsing) : type_scope.
Notation "''M[' R ]_ ( n )" := 'M[R]_n (only parsing) : type_scope.
Notation "''M_' ( m , n )" := 'M[_]_(m, n) : type_scope.
Notation "''rV_' n" := 'M_(1, n) : type_scope.
Notation "''cV_' n" := 'M_(n, 1) : type_scope.
Notation "''M_' n" := 'M_(n, n) : type_scope.
Notation "''M_' ( n )" := 'M_n (only parsing) : type_scope.

Notation "\matrix[ k ]_ ( i , j ) E" := (matrix_of_fun k (fun i jE)) :
   ring_scope.

Notation "\matrix_ ( i < m , j < n ) E" :=
  (@matrix_of_fun _ m n matrix_key (fun i jE)) (only parsing) : ring_scope.

Notation "\matrix_ ( i , j < n ) E" :=
  (\matrix_(i < n, j < n) E) (only parsing) : ring_scope.

Notation "\matrix_ ( i , j ) E" := (\matrix_(i < _, j < _) E) : ring_scope.

Notation "\matrix_ ( i < m ) E" :=
  (\matrix_(i < m, j < _) @fun_of_matrix _ 1 _ E 0 j)
  (only parsing) : ring_scope.
Notation "\matrix_ i E" := (\matrix_(i < _) E) : ring_scope.

Notation "\col_ ( i < n ) E" := (@matrix_of_fun _ n 1 matrix_key (fun i _E))
  (only parsing) : ring_scope.
Notation "\col_ i E" := (\col_(i < _) E) : ring_scope.

Notation "\row_ ( j < n ) E" := (@matrix_of_fun _ 1 n matrix_key (fun _ jE))
  (only parsing) : ring_scope.
Notation "\row_ j E" := (\row_(j < _) E) : ring_scope.


Lemma card_mx (F : finType) m n : (#|{: 'M[F]_(m, n)}| = #|F| ^ (m × n))%N.

#[deprecated(since="mathcomp 1.13.0", note="Use card_mx instead.")]
Notation card_matrix := card_mx.

Matrix structural operations (transpose, permutation, blocks) *******

Section MatrixStructural.

Variable R : Type.

Constant matrix
Fact const_mx_key : unit.
Definition const_mx m n a : 'M[R]_(m, n) := \matrix[const_mx_key]_(i, j) a.
Arguments const_mx {m n}.

Section FixedDim.
Definitions and properties for which we can work with fixed dimensions.

Variables m n : nat.
Implicit Type A : 'M[R]_(m, n).

Reshape a matrix, to accomodate the block functions for instance.
Definition castmx m' n' (eq_mn : (m = m') × (n = n')) A : 'M_(m', n') :=
  let: erefl in _ = m' := eq_mn.1 return 'M_(m', n') in
  let: erefl in _ = n' := eq_mn.2 return 'M_(m, n') in A.

Definition conform_mx m' n' B A :=
  match m =P m', n =P n' with
  | ReflectT eq_m, ReflectT eq_ncastmx (eq_m, eq_n) A
  | _, _B
  end.

Transpose a matrix
Fact trmx_key : unit.
Definition trmx A := \matrix[trmx_key]_(i, j) A j i.

Permute a matrix vertically (rows) or horizontally (columns)
Fact row_perm_key : unit.
Definition row_perm (s : 'S_m) A := \matrix[row_perm_key]_(i, j) A (s i) j.
Fact col_perm_key : unit.
Definition col_perm (s : 'S_n) A := \matrix[col_perm_key]_(i, j) A i (s j).

Exchange two rows/columns of a matrix
Definition xrow i1 i2 := row_perm (tperm i1 i2).
Definition xcol j1 j2 := col_perm (tperm j1 j2).

Row/Column sub matrices of a matrix
Definition row i0 A := \row_j A i0 j.
Definition col j0 A := \col_i A i j0.

Removing a row/column from a matrix
Definition row' i0 A := \matrix_(i, j) A (lift i0 i) j.
Definition col' j0 A := \matrix_(i, j) A i (lift j0 j).

reindexing/subindex a matrix
Definition mxsub m' n' f g A := \matrix_(i < m', j < n') A (f i) (g j).

Lemma castmx_const m' n' (eq_mn : (m = m') × (n = n')) a :
  castmx eq_mn (const_mx a) = const_mx a.

Lemma trmx_const a : trmx (const_mx a) = const_mx a.

Lemma row_perm_const s a : row_perm s (const_mx a) = const_mx a.

Lemma col_perm_const s a : col_perm s (const_mx a) = const_mx a.

Lemma xrow_const i1 i2 a : xrow i1 i2 (const_mx a) = const_mx a.

Lemma xcol_const j1 j2 a : xcol j1 j2 (const_mx a) = const_mx a.

Lemma rowP (u v : 'rV[R]_n) : u 0 =1 v 0 u = v.

Lemma rowK u_ i0 : row i0 (\matrix_i u_ i) = u_ i0.

Lemma row_matrixP A B : ( i, row i A = row i B) A = B.

Lemma colP (u v : 'cV[R]_m) : u^~ 0 =1 v^~ 0 u = v.

Lemma row_const i0 a : row i0 (const_mx a) = const_mx a.

Lemma col_const j0 a : col j0 (const_mx a) = const_mx a.

Lemma row'_const i0 a : row' i0 (const_mx a) = const_mx a.

Lemma col'_const j0 a : col' j0 (const_mx a) = const_mx a.

Lemma col_perm1 A : col_perm 1 A = A.

Lemma row_perm1 A : row_perm 1 A = A.

Lemma col_permM s t A : col_perm (s × t) A = col_perm s (col_perm t A).

Lemma row_permM s t A : row_perm (s × t) A = row_perm s (row_perm t A).

Lemma col_row_permC s t A :
  col_perm s (row_perm t A) = row_perm t (col_perm s A).

Lemma rowEsub i : row i = rowsub (fun i).
Lemma colEsub j : col j = colsub (fun j).

Lemma row'Esub i : row' i = rowsub (lift i).
Lemma col'Esub j : col' j = colsub (lift j).

Lemma row_permEsub s : row_perm s = rowsub s.

Lemma col_permEsub s : col_perm s = colsub s.

Lemma xrowEsub i1 i2 : xrow i1 i2 = rowsub (tperm i1 i2).

Lemma xcolEsub j1 j2 : xcol j1 j2 = colsub (tperm j1 j2).

Lemma mxsub_id : mxsub id id =1 id.

Lemma eq_mxsub m' n' f f' g g' : f =1 f' g =1 g'
  @mxsub m' n' f g =1 mxsub f' g'.

Lemma eq_rowsub m' (f f' : 'I_m' 'I_m) : f =1 f' rowsub f =1 rowsub f'.

Lemma eq_colsub n' (g g' : 'I_n' 'I_n) : g =1 g' colsub g =1 colsub g'.

Lemma mxsub_eq_id f g : f =1 id g =1 id mxsub f g =1 id.

Lemma mxsub_eq_colsub n' f g : f =1 id @mxsub _ n' f g =1 colsub g.

Lemma mxsub_eq_rowsub m' f g : g =1 id @mxsub m' _ f g =1 rowsub f.

Lemma mxsub_ffunl m' n' f g : @mxsub m' n' (finfun f) g =1 mxsub f g.

Lemma mxsub_ffunr m' n' f g : @mxsub m' n' f (finfun g) =1 mxsub f g.

Lemma mxsub_ffun m' n' f g : @mxsub m' n' (finfun f) (finfun g) =1 mxsub f g.

Lemma mxsub_const m' n' f g a : @mxsub m' n' f g (const_mx a) = const_mx a.

End FixedDim.


Lemma castmx_id m n erefl_mn (A : 'M_(m, n)) : castmx erefl_mn A = A.

Lemma castmx_comp m1 n1 m2 n2 m3 n3 (eq_m1 : m1 = m2) (eq_n1 : n1 = n2)
                                    (eq_m2 : m2 = m3) (eq_n2 : n2 = n3) A :
  castmx (eq_m2, eq_n2) (castmx (eq_m1, eq_n1) A)
    = castmx (etrans eq_m1 eq_m2, etrans eq_n1 eq_n2) A.

Lemma castmxK m1 n1 m2 n2 (eq_m : m1 = m2) (eq_n : n1 = n2) :
  cancel (castmx (eq_m, eq_n)) (castmx (esym eq_m, esym eq_n)).

Lemma castmxKV m1 n1 m2 n2 (eq_m : m1 = m2) (eq_n : n1 = n2) :
  cancel (castmx (esym eq_m, esym eq_n)) (castmx (eq_m, eq_n)).

This can be use to reverse an equation that involves a cast.
Lemma castmx_sym m1 n1 m2 n2 (eq_m : m1 = m2) (eq_n : n1 = n2) A1 A2 :
  A1 = castmx (eq_m, eq_n) A2 A2 = castmx (esym eq_m, esym eq_n) A1.

Lemma eq_castmx m1 n1 m2 n2 (eq_mn eq_mn' : (m1 = m2) × (n1 = n2)) :
  castmx eq_mn =1 castmx eq_mn'.

Lemma castmxE m1 n1 m2 n2 (eq_mn : (m1 = m2) × (n1 = n2)) A i j :
  castmx eq_mn A i j =
     A (cast_ord (esym eq_mn.1) i) (cast_ord (esym eq_mn.2) j).

Lemma conform_mx_id m n (B A : 'M_(m, n)) : conform_mx B A = A.

Lemma nonconform_mx m m' n n' (B : 'M_(m', n')) (A : 'M_(m, n)) :
  (m != m') || (n != n') conform_mx B A = B.

Lemma conform_castmx m1 n1 m2 n2 m3 n3
                     (e_mn : (m2 = m3) × (n2 = n3)) (B : 'M_(m1, n1)) A :
  conform_mx B (castmx e_mn A) = conform_mx B A.

Lemma trmxK m n : cancel (@trmx m n) (@trmx n m).

Lemma trmx_inj m n : injective (@trmx m n).

Lemma trmx_cast m1 n1 m2 n2 (eq_mn : (m1 = m2) × (n1 = n2)) A :
  (castmx eq_mn A)^T = castmx (eq_mn.2, eq_mn.1) A^T.

Lemma trmx_conform m' n' m n (B : 'M_(m', n')) (A : 'M_(m, n)) :
  (conform_mx B A)^T = conform_mx B^T A^T.

Lemma tr_row_perm m n s (A : 'M_(m, n)) : (row_perm s A)^T = col_perm s A^T.

Lemma tr_col_perm m n s (A : 'M_(m, n)) : (col_perm s A)^T = row_perm s A^T.

Lemma tr_xrow m n i1 i2 (A : 'M_(m, n)) : (xrow i1 i2 A)^T = xcol i1 i2 A^T.

Lemma tr_xcol m n j1 j2 (A : 'M_(m, n)) : (xcol j1 j2 A)^T = xrow j1 j2 A^T.

Lemma row_id n i (V : 'rV_n) : row i V = V.

Lemma col_id n j (V : 'cV_n) : col j V = V.

Lemma row_eq m1 m2 n i1 i2 (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) :
  row i1 A1 = row i2 A2 A1 i1 =1 A2 i2.

Lemma col_eq m n1 n2 j1 j2 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :
  col j1 A1 = col j2 A2 A1^~ j1 =1 A2^~ j2.

Lemma row'_eq m n i0 (A B : 'M_(m, n)) :
  row' i0 A = row' i0 B {in predC1 i0, A =2 B}.

Lemma col'_eq m n j0 (A B : 'M_(m, n)) :
  col' j0 A = col' j0 B i, {in predC1 j0, A i =1 B i}.

Lemma tr_row m n i0 (A : 'M_(m, n)) : (row i0 A)^T = col i0 A^T.

Lemma tr_row' m n i0 (A : 'M_(m, n)) : (row' i0 A)^T = col' i0 A^T.

Lemma tr_col m n j0 (A : 'M_(m, n)) : (col j0 A)^T = row j0 A^T.

Lemma tr_col' m n j0 (A : 'M_(m, n)) : (col' j0 A)^T = row' j0 A^T.

Lemma mxsub_comp m1 m2 m3 n1 n2 n3
  (f : 'I_m2 'I_m1) (f' : 'I_m3 'I_m2)
  (g : 'I_n2 'I_n1) (g' : 'I_n3 'I_n2) (A : 'M_(m1, n1)) :
  mxsub (f \o f') (g \o g') A = mxsub f' g' (mxsub f g A).

Lemma rowsub_comp m1 m2 m3 n
  (f : 'I_m2 'I_m1) (f' : 'I_m3 'I_m2) (A : 'M_(m1, n)) :
  rowsub (f \o f') A = rowsub f' (rowsub f A).

Lemma colsub_comp m n n2 n3
  (g : 'I_n2 'I_n) (g' : 'I_n3 'I_n2) (A : 'M_(m, n)) :
  colsub (g \o g') A = colsub g' (colsub g A).

Lemma mxsubrc m1 m2 n n2 f g (A : 'M_(m1, n)) :
  mxsub f g A = rowsub f (colsub g A) :> 'M_(m2, n2).

Lemma mxsubcr m1 m2 n n2 f g (A : 'M_(m1, n)) :
  mxsub f g A = colsub g (rowsub f A) :> 'M_(m2, n2).

Lemma rowsub_cast m1 m2 n (eq_m : m1 = m2) (A : 'M_(m2, n)) :
  rowsub (cast_ord eq_m) A = castmx (esym eq_m, erefl) A.

Lemma colsub_cast m n1 n2 (eq_n : n1 = n2) (A : 'M_(m, n2)) :
  colsub (cast_ord eq_n) A = castmx (erefl, esym eq_n) A.

Lemma mxsub_cast m1 m2 n1 n2 (eq_m : m1 = m2) (eq_n : n1 = n2) A :
  mxsub (cast_ord eq_m) (cast_ord eq_n) A = castmx (esym eq_m, esym eq_n) A.

Lemma castmxEsub m1 m2 n1 n2 (eq_mn : (m1 = m2) × (n1 = n2)) A :
  castmx eq_mn A = mxsub (cast_ord (esym eq_mn.1)) (cast_ord (esym eq_mn.2)) A.

Lemma trmx_mxsub m1 m2 n1 n2 f g (A : 'M_(m1, n1)) :
  (mxsub f g A)^T = mxsub g f A^T :> 'M_(n2, m2).

Lemma row_mxsub m1 m2 n1 n2
    (f : 'I_m2 'I_m1) (g : 'I_n2 'I_n1) (A : 'M_(m1, n1)) i :
  row i (mxsub f g A) = row (f i) (colsub g A).

Lemma col_mxsub m1 m2 n1 n2
    (f : 'I_m2 'I_m1) (g : 'I_n2 'I_n1) (A : 'M_(m1, n1)) i :
 col i (mxsub f g A) = col (g i) (rowsub f A).

Lemma row_rowsub m1 m2 n (f : 'I_m2 'I_m1) (A : 'M_(m1, n)) i :
  row i (rowsub f A) = row (f i) A.

Lemma col_colsub m n1 n2 (g : 'I_n2 'I_n1) (A : 'M_(m, n1)) i :
  col i (colsub g A) = col (g i) A.

Ltac split_mxE := apply/matrixPi j; do ![rewrite mxE | case: split ⇒ ?].

Section CutPaste.

Variables m m1 m2 n n1 n2 : nat.

Concatenating two matrices, in either direction.

Fact row_mx_key : unit.
Definition row_mx (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) : 'M[R]_(m, n1 + n2) :=
  \matrix[row_mx_key]_(i, j)
     match split j with inl j1A1 i j1 | inr j2A2 i j2 end.

Fact col_mx_key : unit.
Definition col_mx (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) : 'M[R]_(m1 + m2, n) :=
  \matrix[col_mx_key]_(i, j)
     match split i with inl i1A1 i1 j | inr i2A2 i2 j end.

Left/Right | Up/Down submatrices of a rows | columns matrix. The shape of the (dependent) width parameters of the type of A determines which submatrix is selected.

Fact lsubmx_key : unit.
Definition lsubmx (A : 'M[R]_(m, n1 + n2)) :=
  \matrix[lsubmx_key]_(i, j) A i (lshift n2 j).

Fact rsubmx_key : unit.
Definition rsubmx (A : 'M[R]_(m, n1 + n2)) :=
  \matrix[rsubmx_key]_(i, j) A i (rshift n1 j).

Fact usubmx_key : unit.
Definition usubmx (A : 'M[R]_(m1 + m2, n)) :=
  \matrix[usubmx_key]_(i, j) A (lshift m2 i) j.

Fact dsubmx_key : unit.
Definition dsubmx (A : 'M[R]_(m1 + m2, n)) :=
  \matrix[dsubmx_key]_(i, j) A (rshift m1 i) j.

Lemma row_mxEl A1 A2 i j : row_mx A1 A2 i (lshift n2 j) = A1 i j.

Lemma row_mxKl A1 A2 : lsubmx (row_mx A1 A2) = A1.

Lemma row_mxEr A1 A2 i j : row_mx A1 A2 i (rshift n1 j) = A2 i j.

Lemma row_mxKr A1 A2 : rsubmx (row_mx A1 A2) = A2.

Lemma hsubmxK A : row_mx (lsubmx A) (rsubmx A) = A.

Lemma col_mxEu A1 A2 i j : col_mx A1 A2 (lshift m2 i) j = A1 i j.

Lemma col_mxKu A1 A2 : usubmx (col_mx A1 A2) = A1.

Lemma col_mxEd A1 A2 i j : col_mx A1 A2 (rshift m1 i) j = A2 i j.

Lemma col_mxKd A1 A2 : dsubmx (col_mx A1 A2) = A2.

Lemma lsubmxEsub : lsubmx = colsub (lshift _).

Lemma rsubmxEsub : rsubmx = colsub (@rshift _ _).

Lemma usubmxEsub : usubmx = rowsub (lshift _).

Lemma dsubmxEsub : dsubmx = rowsub (@rshift _ _).

Lemma eq_row_mx A1 A2 B1 B2 : row_mx A1 A2 = row_mx B1 B2 A1 = B1 A2 = B2.

Lemma eq_col_mx A1 A2 B1 B2 : col_mx A1 A2 = col_mx B1 B2 A1 = B1 A2 = B2.

Lemma row_mx_const a : row_mx (const_mx a) (const_mx a) = const_mx a.

Lemma col_mx_const a : col_mx (const_mx a) (const_mx a) = const_mx a.

Lemma row_usubmx A i : row i (usubmx A) = row (lshift m2 i) A.

Lemma row_dsubmx A i : row i (dsubmx A) = row (rshift m1 i) A.

Lemma col_lsubmx A i : col i (lsubmx A) = col (lshift n2 i) A.

Lemma col_rsubmx A i : col i (rsubmx A) = col (rshift n1 i) A.

End CutPaste.

Lemma row_thin_mx m n (A : 'M_(m,0)) (B : 'M_(m,n)) : row_mx A B = B.

Lemma col_flat_mx m n (A : 'M_(0,n)) (B : 'M_(m,n)) : col_mx A B = B.

Lemma trmx_lsub m n1 n2 (A : 'M_(m, n1 + n2)) : (lsubmx A)^T = usubmx A^T.

Lemma trmx_rsub m n1 n2 (A : 'M_(m, n1 + n2)) : (rsubmx A)^T = dsubmx A^T.

Lemma tr_row_mx m n1 n2 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :
  (row_mx A1 A2)^T = col_mx A1^T A2^T.

Lemma tr_col_mx m1 m2 n (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) :
  (col_mx A1 A2)^T = row_mx A1^T A2^T.

Lemma trmx_usub m1 m2 n (A : 'M_(m1 + m2, n)) : (usubmx A)^T = lsubmx A^T.

Lemma trmx_dsub m1 m2 n (A : 'M_(m1 + m2, n)) : (dsubmx A)^T = rsubmx A^T.

Lemma vsubmxK m1 m2 n (A : 'M_(m1 + m2, n)) : col_mx (usubmx A) (dsubmx A) = A.

Lemma cast_row_mx m m' n1 n2 (eq_m : m = m') A1 A2 :
  castmx (eq_m, erefl _) (row_mx A1 A2)
    = row_mx (castmx (eq_m, erefl n1) A1) (castmx (eq_m, erefl n2) A2).

Lemma cast_col_mx m1 m2 n n' (eq_n : n = n') A1 A2 :
  castmx (erefl _, eq_n) (col_mx A1 A2)
    = col_mx (castmx (erefl m1, eq_n) A1) (castmx (erefl m2, eq_n) A2).

This lemma has Prenex Implicits to help RL rewrititng with castmx_sym.
Lemma row_mxA m n1 n2 n3 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) (A3 : 'M_(m, n3)) :
  let cast := (erefl m, esym (addnA n1 n2 n3)) in
  row_mx A1 (row_mx A2 A3) = castmx cast (row_mx (row_mx A1 A2) A3).
Definition row_mxAx := row_mxA. (* bypass Prenex Implicits. *)

This lemma has Prenex Implicits to help RL rewrititng with castmx_sym.
Lemma col_mxA m1 m2 m3 n (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) (A3 : 'M_(m3, n)) :
  let cast := (esym (addnA m1 m2 m3), erefl n) in
  col_mx A1 (col_mx A2 A3) = castmx cast (col_mx (col_mx A1 A2) A3).
Definition col_mxAx := col_mxA. (* bypass Prenex Implicits. *)

Lemma row_row_mx m n1 n2 i0 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :
  row i0 (row_mx A1 A2) = row_mx (row i0 A1) (row i0 A2).

Lemma col_col_mx m1 m2 n j0 (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) :
  col j0 (col_mx A1 A2) = col_mx (col j0 A1) (col j0 A2).

Lemma row'_row_mx m n1 n2 i0 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :
  row' i0 (row_mx A1 A2) = row_mx (row' i0 A1) (row' i0 A2).

Lemma col'_col_mx m1 m2 n j0 (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) :
  col' j0 (col_mx A1 A2) = col_mx (col' j0 A1) (col' j0 A2).

Lemma colKl m n1 n2 j1 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :
  col (lshift n2 j1) (row_mx A1 A2) = col j1 A1.

Lemma colKr m n1 n2 j2 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :
  col (rshift n1 j2) (row_mx A1 A2) = col j2 A2.

Lemma rowKu m1 m2 n i1 (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) :
  row (lshift m2 i1) (col_mx A1 A2) = row i1 A1.

Lemma rowKd m1 m2 n i2 (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) :
  row (rshift m1 i2) (col_mx A1 A2) = row i2 A2.

Lemma col'Kl m n1 n2 j1 (A1 : 'M_(m, n1.+1)) (A2 : 'M_(m, n2)) :
  col' (lshift n2 j1) (row_mx A1 A2) = row_mx (col' j1 A1) A2.

Lemma row'Ku m1 m2 n i1 (A1 : 'M_(m1.+1, n)) (A2 : 'M_(m2, n)) :
  row' (lshift m2 i1) (@col_mx m1.+1 m2 n A1 A2) = col_mx (row' i1 A1) A2.

Lemma mx'_cast m n : 'I_n (m + n.-1)%N = (m + n).-1.

Lemma col'Kr m n1 n2 j2 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :
  col' (rshift n1 j2) (@row_mx m n1 n2 A1 A2)
    = castmx (erefl m, mx'_cast n1 j2) (row_mx A1 (col' j2 A2)).

Lemma row'Kd m1 m2 n i2 (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) :
  row' (rshift m1 i2) (col_mx A1 A2)
    = castmx (mx'_cast m1 i2, erefl n) (col_mx A1 (row' i2 A2)).

Section Block.

Variables m1 m2 n1 n2 : nat.

Building a block matrix from 4 matrices : up left, up right, down left and down right components

Definition block_mx Aul Aur Adl Adr : 'M_(m1 + m2, n1 + n2) :=
  col_mx (row_mx Aul Aur) (row_mx Adl Adr).

Lemma eq_block_mx Aul Aur Adl Adr Bul Bur Bdl Bdr :
 block_mx Aul Aur Adl Adr = block_mx Bul Bur Bdl Bdr
  [/\ Aul = Bul, Aur = Bur, Adl = Bdl & Adr = Bdr].

Lemma block_mx_const a :
  block_mx (const_mx a) (const_mx a) (const_mx a) (const_mx a) = const_mx a.

Section CutBlock.

Variable A : matrix R (m1 + m2) (n1 + n2).

Definition ulsubmx := lsubmx (usubmx A).
Definition ursubmx := rsubmx (usubmx A).
Definition dlsubmx := lsubmx (dsubmx A).
Definition drsubmx := rsubmx (dsubmx A).

Lemma submxK : block_mx ulsubmx ursubmx dlsubmx drsubmx = A.

Lemma ulsubmxEsub : ulsubmx = mxsub (lshift _) (lshift _) A.

Lemma dlsubmxEsub : dlsubmx = mxsub (@rshift _ _) (lshift _) A.

Lemma ursubmxEsub : ursubmx = mxsub (lshift _) (@rshift _ _) A.

Lemma drsubmxEsub : drsubmx = mxsub (@rshift _ _) (@rshift _ _) A.

End CutBlock.

Section CatBlock.

Variables (Aul : 'M[R]_(m1, n1)) (Aur : 'M[R]_(m1, n2)).
Variables (Adl : 'M[R]_(m2, n1)) (Adr : 'M[R]_(m2, n2)).

Let A := block_mx Aul Aur Adl Adr.

Lemma block_mxEul i j : A (lshift m2 i) (lshift n2 j) = Aul i j.
Lemma block_mxKul : ulsubmx A = Aul.

Lemma block_mxEur i j : A (lshift m2 i) (rshift n1 j) = Aur i j.
Lemma block_mxKur : ursubmx A = Aur.

Lemma block_mxEdl i j : A (rshift m1 i) (lshift n2 j) = Adl i j.
Lemma block_mxKdl : dlsubmx A = Adl.

Lemma block_mxEdr i j : A (rshift m1 i) (rshift n1 j) = Adr i j.
Lemma block_mxKdr : drsubmx A = Adr.

Lemma block_mxEv : A = col_mx (row_mx Aul Aur) (row_mx Adl Adr).

End CatBlock.

End Block.

Section TrCutBlock.

Variables m1 m2 n1 n2 : nat.
Variable A : 'M[R]_(m1 + m2, n1 + n2).

Lemma trmx_ulsub : (ulsubmx A)^T = ulsubmx A^T.

Lemma trmx_ursub : (ursubmx A)^T = dlsubmx A^T.

Lemma trmx_dlsub : (dlsubmx A)^T = ursubmx A^T.

Lemma trmx_drsub : (drsubmx A)^T = drsubmx A^T.

End TrCutBlock.

Section TrBlock.
Variables m1 m2 n1 n2 : nat.
Variables (Aul : 'M[R]_(m1, n1)) (Aur : 'M[R]_(m1, n2)).
Variables (Adl : 'M[R]_(m2, n1)) (Adr : 'M[R]_(m2, n2)).

Lemma tr_block_mx :
 (block_mx Aul Aur Adl Adr)^T = block_mx Aul^T Adl^T Aur^T Adr^T.

Lemma block_mxEh :
  block_mx Aul Aur Adl Adr = row_mx (col_mx Aul Adl) (col_mx Aur Adr).
End TrBlock.

This lemma has Prenex Implicits to help RL rewrititng with castmx_sym.
Lemma block_mxA m1 m2 m3 n1 n2 n3
   (A11 : 'M_(m1, n1)) (A12 : 'M_(m1, n2)) (A13 : 'M_(m1, n3))
   (A21 : 'M_(m2, n1)) (A22 : 'M_(m2, n2)) (A23 : 'M_(m2, n3))
   (A31 : 'M_(m3, n1)) (A32 : 'M_(m3, n2)) (A33 : 'M_(m3, n3)) :
  let cast := (esym (addnA m1 m2 m3), esym (addnA n1 n2 n3)) in
  let row1 := row_mx A12 A13 in let col1 := col_mx A21 A31 in
  let row3 := row_mx A31 A32 in let col3 := col_mx A13 A23 in
  block_mx A11 row1 col1 (block_mx A22 A23 A32 A33)
    = castmx cast (block_mx (block_mx A11 A12 A21 A22) col3 row3 A33).
Definition block_mxAx := block_mxA. (* Bypass Prenex Implicits *)

Section Induction.

Lemma row_ind m (P : n, 'M[R]_(m, n) Type) :
    ( A, P 0%N A)
    ( n c A, P n A P (1 + n)%N (row_mx c A))
   n A, P n A.

Lemma col_ind n (P : m, 'M[R]_(m, n) Type) :
    ( A, P 0%N A)
    ( m r A, P m A P (1 + m)%N (col_mx r A))
   m A, P m A.

Lemma mx_ind (P : m n, 'M[R]_(m, n) Type) :
    ( m A, P m 0%N A)
    ( n A, P 0%N n A)
    ( m n x r c A, P m n A P (1 + m)%N (1 + n)%N (block_mx x r c A))
   m n A, P m n A.
Definition matrix_rect := mx_ind.
Definition matrix_rec := mx_ind.
Definition matrix_ind := mx_ind.

Lemma sqmx_ind (P : n, 'M[R]_n Type) :
    ( A, P 0%N A)
    ( n x r c A, P n A P (1 + n)%N (block_mx x r c A))
   n A, P n A.

Lemma ringmx_ind (P : n, 'M[R]_n.+1 Type) :
    ( x, P 0%N x)
    ( n x (r : 'rV_n.+1) (c : 'cV_n.+1) A,
       P n A P (1 + n)%N (block_mx x r c A))
   n A, P n A.

Lemma mxsub_ind
    (weight : m n, 'M[R]_(m, n) nat)
    (sub : m n m' n', ('I_m' 'I_m) ('I_n' 'I_n) Prop)
    (P : m n, 'M[R]_(m, n) Type) :
    ( m n (A : 'M[R]_(m, n)),
      ( m' n' f g, weight m' n' (mxsub f g A) < weight m n A
                         sub m n m' n' f g
                         P m' n' (mxsub f g A)) P m n A)
   m n A, P m n A.

End Induction.

Bijections mxvec : 'M(m, n) <----> 'rV(m * n) : vec_mx
Section VecMatrix.

Variables m n : nat.

Lemma mxvec_cast : #|{:'I_m × 'I_n}| = (m × n)%N.

Definition mxvec_index (i : 'I_m) (j : 'I_n) :=
  cast_ord mxvec_cast (enum_rank (i, j)).

Variant is_mxvec_index : 'I_(m × n) Type :=
  isMxvecIndex i j : is_mxvec_index (mxvec_index i j).

Lemma mxvec_indexP k : is_mxvec_index k.

Coercion pair_of_mxvec_index k (i_k : is_mxvec_index k) :=
  let: isMxvecIndex i j := i_k in (i, j).

Definition mxvec (A : 'M[R]_(m, n)) :=
  castmx (erefl _, mxvec_cast) (\row_k A (enum_val k).1 (enum_val k).2).

Fact vec_mx_key : unit.
Definition vec_mx (u : 'rV[R]_(m × n)) :=
  \matrix[vec_mx_key]_(i, j) u 0 (mxvec_index i j).

Lemma mxvecE A i j : mxvec A 0 (mxvec_index i j) = A i j.

Lemma mxvecK : cancel mxvec vec_mx.

Lemma vec_mxK : cancel vec_mx mxvec.

Lemma curry_mxvec_bij : {on 'I_(m × n), bijective (uncurry mxvec_index)}.

End VecMatrix.

End MatrixStructural.

Arguments const_mx {R m n}.
Arguments row_mxA {R m n1 n2 n3 A1 A2 A3}.
Arguments col_mxA {R m1 m2 m3 n A1 A2 A3}.
Arguments block_mxA
  {R m1 m2 m3 n1 n2 n3 A11 A12 A13 A21 A22 A23 A31 A32 A33}.
Arguments trmx_inj {R m n} [A1 A2] eqA12t : rename.

Notation "A ^T" := (trmx A) : ring_scope.
Notation colsub g := (mxsub id g).
Notation rowsub f := (mxsub f id).

Arguments eq_mxsub [R m n m' n' f] f' [g] g' _.
Arguments eq_rowsub [R m n m' f] f' _.
Arguments eq_colsub [R m n n' g] g' _.

Matrix parametricity.
Section MapMatrix.

Variables (aT rT : Type) (f : aT rT).

Fact map_mx_key : unit.
Definition map_mx m n (A : 'M_(m, n)) := \matrix[map_mx_key]_(i, j) f (A i j).

Notation "A ^f" := (map_mx A) : ring_scope.

Section OneMatrix.

Variables (m n : nat) (A : 'M[aT]_(m, n)).

Lemma map_trmx : A^f^T = A^T^f.

Lemma map_const_mx a : (const_mx a)^f = const_mx (f a) :> 'M_(m, n).

Lemma map_row i : (row i A)^f = row i A^f.

Lemma map_col j : (col j A)^f = col j A^f.

Lemma map_row' i0 : (row' i0 A)^f = row' i0 A^f.

Lemma map_col' j0 : (col' j0 A)^f = col' j0 A^f.

Lemma map_mxsub m' n' g h : (@mxsub _ _ _ m' n' g h A)^f = mxsub g h A^f.

Lemma map_row_perm s : (row_perm s A)^f = row_perm s A^f.

Lemma map_col_perm s : (col_perm s A)^f = col_perm s A^f.

Lemma map_xrow i1 i2 : (xrow i1 i2 A)^f = xrow i1 i2 A^f.

Lemma map_xcol j1 j2 : (xcol j1 j2 A)^f = xcol j1 j2 A^f.

Lemma map_castmx m' n' c : (castmx c A)^f = castmx c A^f :> 'M_(m', n').

Lemma map_conform_mx m' n' (B : 'M_(m', n')) :
  (conform_mx B A)^f = conform_mx B^f A^f.

Lemma map_mxvec : (mxvec A)^f = mxvec A^f.

Lemma map_vec_mx (v : 'rV_(m × n)) : (vec_mx v)^f = vec_mx v^f.

End OneMatrix.

Section Block.

Variables m1 m2 n1 n2 : nat.
Variables (Aul : 'M[aT]_(m1, n1)) (Aur : 'M[aT]_(m1, n2)).
Variables (Adl : 'M[aT]_(m2, n1)) (Adr : 'M[aT]_(m2, n2)).
Variables (Bh : 'M[aT]_(m1, n1 + n2)) (Bv : 'M[aT]_(m1 + m2, n1)).
Variable B : 'M[aT]_(m1 + m2, n1 + n2).

Lemma map_row_mx : (row_mx Aul Aur)^f = row_mx Aul^f Aur^f.

Lemma map_col_mx : (col_mx Aul Adl)^f = col_mx Aul^f Adl^f.

Lemma map_block_mx :
  (block_mx Aul Aur Adl Adr)^f = block_mx Aul^f Aur^f Adl^f Adr^f.

Lemma map_lsubmx : (lsubmx Bh)^f = lsubmx Bh^f.

Lemma map_rsubmx : (rsubmx Bh)^f = rsubmx Bh^f.

Lemma map_usubmx : (usubmx Bv)^f = usubmx Bv^f.

Lemma map_dsubmx : (dsubmx Bv)^f = dsubmx Bv^f.

Lemma map_ulsubmx : (ulsubmx B)^f = ulsubmx B^f.

Lemma map_ursubmx : (ursubmx B)^f = ursubmx B^f.

Lemma map_dlsubmx : (dlsubmx B)^f = dlsubmx B^f.

Lemma map_drsubmx : (drsubmx B)^f = drsubmx B^f.

End Block.

End MapMatrix.

Arguments map_mx {aT rT} f {m n} A.

Section MultipleMapMatrix.
Context {R S T : Type} {m n : nat}.

Lemma map_mx_comp (f : R S) (g : S T)
  (M : 'M_(m, n)) : M ^ (g \o f) = (M ^ f) ^ g.

Lemma eq_in_map_mx (g f : R S) (M : 'M_(m, n)) :
  ( i j, f (M i j) = g (M i j)) M ^ f = M ^ g.

Lemma eq_map_mx (g f : R S) : f =1 g
   (M : 'M_(m, n)), M ^ f = M ^ g.

Lemma map_mx_id_in (f : R R) (M : 'M_(m, n)) :
  ( i j, f (M i j) = M i j) M ^ f = M.

Lemma map_mx_id (f : R R) : f =1 id M : 'M_(m, n), M ^ f = M.

End MultipleMapMatrix.
Arguments eq_map_mx {R S m n} g [f].
Arguments eq_in_map_mx {R S m n} g [f M].
Arguments map_mx_id_in {R m n} [f M].
Arguments map_mx_id {R m n} [f].

Matrix Zmodule (additive) structure ******************

Section MatrixZmodule.

Variable V : zmodType.

Section FixedDim.

Variables m n : nat.
Implicit Types A B : 'M[V]_(m, n).

Fact oppmx_key : unit.
Fact addmx_key : unit.
Definition oppmx A := \matrix[oppmx_key]_(i, j) (- A i j).
Definition addmx A B := \matrix[addmx_key]_(i, j) (A i j + B i j).
In principle, diag_mx and scalar_mx could be defined here, but since they only make sense with the graded ring operations, we defer them to the next section.

Lemma addmxA : associative addmx.

Lemma addmxC : commutative addmx.

Lemma add0mx : left_id (const_mx 0) addmx.

Lemma addNmx : left_inverse (const_mx 0) oppmx addmx.


Lemma mulmxnE A d i j : (A *+ d) i j = A i j *+ d.

Lemma summxE I r (P : pred I) (E : I 'M_(m, n)) i j :
  (\sum_(k <- r | P k) E k) i j = \sum_(k <- r | P k) E k i j.

Lemma const_mx_is_additive : additive const_mx.

End FixedDim.

Section Additive.

Variables (m n p q : nat) (f : 'I_p 'I_q 'I_m) (g : 'I_p 'I_q 'I_n).

Definition swizzle_mx k (A : 'M[V]_(m, n)) :=
  \matrix[k]_(i, j) A (f i j) (g i j).

Lemma swizzle_mx_is_additive k : additive (swizzle_mx k).

End Additive.



Lemma flatmx0 n : all_equal_to (0 : 'M_(0, n)).

Lemma thinmx0 n : all_equal_to (0 : 'M_(n, 0)).

Lemma trmx0 m n : (0 : 'M_(m, n))^T = 0.

Lemma row0 m n i0 : row i0 (0 : 'M_(m, n)) = 0.

Lemma col0 m n j0 : col j0 (0 : 'M_(m, n)) = 0.

Lemma mxvec_eq0 m n (A : 'M_(m, n)) : (mxvec A == 0) = (A == 0).

Lemma vec_mx_eq0 m n (v : 'rV_(m × n)) : (vec_mx v == 0) = (v == 0).

Lemma row_mx0 m n1 n2 : row_mx 0 0 = 0 :> 'M_(m, n1 + n2).

Lemma col_mx0 m1 m2 n : col_mx 0 0 = 0 :> 'M_(m1 + m2, n).

Lemma block_mx0 m1 m2 n1 n2 : block_mx 0 0 0 0 = 0 :> 'M_(m1 + m2, n1 + n2).

Ltac split_mxE := apply/matrixPi j; do ![rewrite mxE | case: split ⇒ ?].

Lemma opp_row_mx m n1 n2 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :
  - row_mx A1 A2 = row_mx (- A1) (- A2).

Lemma opp_col_mx m1 m2 n (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) :
  - col_mx A1 A2 = col_mx (- A1) (- A2).

Lemma opp_block_mx m1 m2 n1 n2 (Aul : 'M_(m1, n1)) Aur Adl (Adr : 'M_(m2, n2)) :
  - block_mx Aul Aur Adl Adr = block_mx (- Aul) (- Aur) (- Adl) (- Adr).

Lemma add_row_mx m n1 n2 (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) B1 B2 :
  row_mx A1 A2 + row_mx B1 B2 = row_mx (A1 + B1) (A2 + B2).

Lemma add_col_mx m1 m2 n (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) B1 B2 :
  col_mx A1 A2 + col_mx B1 B2 = col_mx (A1 + B1) (A2 + B2).

Lemma add_block_mx m1 m2 n1 n2 (Aul : 'M_(m1, n1)) Aur Adl (Adr : 'M_(m2, n2))
                   Bul Bur Bdl Bdr :
  let A := block_mx Aul Aur Adl Adr in let B := block_mx Bul Bur Bdl Bdr in
  A + B = block_mx (Aul + Bul) (Aur + Bur) (Adl + Bdl) (Adr + Bdr).

Lemma row_mx_eq0 (m n1 n2 : nat) (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)):
  (row_mx A1 A2 == 0) = (A1 == 0) && (A2 == 0).

Lemma col_mx_eq0 (m1 m2 n : nat) (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)):
  (col_mx A1 A2 == 0) = (A1 == 0) && (A2 == 0).

Lemma block_mx_eq0 m1 m2 n1 n2 (Aul : 'M_(m1, n1)) Aur Adl (Adr : 'M_(m2, n2)) :
  (block_mx Aul Aur Adl Adr == 0) =
  [&& Aul == 0, Aur == 0, Adl == 0 & Adr == 0].

Lemma trmx_eq0 m n (A : 'M_(m, n)) : (A^T == 0) = (A == 0).

Lemma matrix_eq0 m n (A : 'M_(m, n)) :
  (A == 0) = [ i, j, A i j == 0].

Lemma matrix0Pn m n (A : 'M_(m, n)) : reflect ( i j, A i j != 0) (A != 0).

Lemma rV0Pn n (v : 'rV_n) : reflect ( i, v 0 i != 0) (v != 0).

Lemma cV0Pn n (v : 'cV_n) : reflect ( i, v i 0 != 0) (v != 0).

Definition nz_row m n (A : 'M_(m, n)) :=
  oapp (fun irow i A) 0 [pick i | row i A != 0].

Lemma nz_row_eq0 m n (A : 'M_(m, n)) : (nz_row A == 0) = (A == 0).

Definition is_diag_mx m n (A : 'M[V]_(m, n)) :=
  [ i : 'I__, j : 'I__, (i != j :> nat) ==> (A i j == 0)].

Lemma is_diag_mxP m n (A : 'M[V]_(m, n)) :
  reflect ( i j : 'I__, i != j :> nat A i j = 0) (is_diag_mx A).

Lemma mx0_is_diag m n : is_diag_mx (0 : 'M[V]_(m, n)).

Lemma mx11_is_diag (M : 'M_1) : is_diag_mx M.

Definition is_trig_mx m n (A : 'M[V]_(m, n)) :=
  [ i : 'I__, j : 'I__, (i < j)%N ==> (A i j == 0)].

Lemma is_trig_mxP m n (A : 'M[V]_(m, n)) :
  reflect ( i j : 'I__, (i < j)%N A i j = 0) (is_trig_mx A).

Lemma is_diag_mx_is_trig m n (A : 'M[V]_(m, n)) : is_diag_mx A is_trig_mx A.

Lemma mx0_is_trig m n : is_trig_mx (0 : 'M[V]_(m, n)).

Lemma mx11_is_trig (M : 'M_1) : is_trig_mx M.

Lemma is_diag_mxEtrig m n (A : 'M[V]_(m, n)) :
  is_diag_mx A = is_trig_mx A && is_trig_mx A^T.

Lemma is_diag_trmx m n (A : 'M[V]_(m, n)) : is_diag_mx A^T = is_diag_mx A.

Lemma ursubmx_trig m1 m2 n1 n2 (A : 'M[V]_(m1 + m2, n1 + n2)) :
  m1 n1 is_trig_mx A ursubmx A = 0.

Lemma dlsubmx_diag m1 m2 n1 n2 (A : 'M[V]_(m1 + m2, n1 + n2)) :
  n1 m1 is_diag_mx A dlsubmx A = 0.

Lemma ulsubmx_trig m1 m2 n1 n2 (A : 'M[V]_(m1 + m2, n1 + n2)) :
  is_trig_mx A is_trig_mx (ulsubmx A).

Lemma drsubmx_trig m1 m2 n1 n2 (A : 'M[V]_(m1 + m2, n1 + n2)) :
  m1 n1 is_trig_mx A is_trig_mx (drsubmx A).

Lemma ulsubmx_diag m1 m2 n1 n2 (A : 'M[V]_(m1 + m2, n1 + n2)) :
  is_diag_mx A is_diag_mx (ulsubmx A).

Lemma drsubmx_diag m1 m2 n1 n2 (A : 'M[V]_(m1 + m2, n1 + n2)) :
  m1 = n1 is_diag_mx A is_diag_mx (drsubmx A).

Lemma is_trig_block_mx m1 m2 n1 n2 ul ur dl dr : m1 = n1
  @is_trig_mx (m1 + m2) (n1 + n2) (block_mx ul ur dl dr) =
  [&& ur == 0, is_trig_mx ul & is_trig_mx dr].

Lemma trigmx_ind (P : m n, 'M_(m, n) Type) :
  ( m, P m 0%N 0)
  ( n, P 0%N n 0)
  ( m n x c A, is_trig_mx A
    P m n A P (1 + m)%N (1 + n)%N (block_mx x 0 c A))
   m n A, is_trig_mx A P m n A.

Lemma trigsqmx_ind (P : n, 'M[V]_n Type) : (P 0%N 0)
  ( n x c A, is_trig_mx A P n A P (1 + n)%N (block_mx x 0 c A))
   n A, is_trig_mx A P n A.

Lemma is_diag_block_mx m1 m2 n1 n2 ul ur dl dr : m1 = n1
  @is_diag_mx (m1 + m2) (n1 + n2) (block_mx ul ur dl dr) =
  [&& ur == 0, dl == 0, is_diag_mx ul & is_diag_mx dr].

Lemma diagmx_ind (P : m n, 'M_(m, n) Type) :
  ( m, P m 0%N 0)
  ( n, P 0%N n 0)
  ( m n x c A, is_diag_mx A
    P m n A P (1 + m)%N (1 + n)%N (block_mx x 0 c A))
   m n A, is_diag_mx A P m n A.

Lemma diagsqmx_ind (P : n, 'M[V]_n Type) :
    (P 0%N 0)
  ( n x c A, is_diag_mx A P n A P (1 + n)%N (block_mx x 0 c A))
   n A, is_diag_mx A P n A.

End MatrixZmodule.

Arguments is_diag_mx {V m n}.
Arguments is_diag_mxP {V m n A}.
Arguments is_trig_mx {V m n}.
Arguments is_trig_mxP {V m n A}.

Parametricity over the additive structure.
Section MapZmodMatrix.

Variables (aR rR : zmodType) (f : {additive aR rR}) (m n : nat).
Implicit Type A : 'M[aR]_(m, n).

Lemma map_mx0 : 0^f = 0 :> 'M_(m, n).

Lemma map_mxN A : (- A)^f = - A^f.

Lemma map_mxD A B : (A + B)^f = A^f + B^f.

Lemma map_mxB A B : (A - B)^f = A^f - B^f.

Definition map_mx_sum := big_morph _ map_mxD map_mx0.


End MapZmodMatrix.

Matrix ring module, graded ring, and ring structures ***********

Section MatrixAlgebra.

Variable R : ringType.

Section RingModule.

The ring module/vector space structure

Variables m n : nat.
Implicit Types A B : 'M[R]_(m, n).

Fact scalemx_key : unit.
Definition scalemx x A := \matrix[scalemx_key]_(i, j) (x × A i j).

Basis
Fact delta_mx_key : unit.
Definition delta_mx i0 j0 : 'M[R]_(m, n) :=
  \matrix[delta_mx_key]_(i, j) ((i == i0) && (j == j0))%:R.


Lemma scale1mx A : 1 ×m: A = A.

Lemma scalemxDl A x y : (x + y) ×m: A = x ×m: A + y ×m: A.

Lemma scalemxDr x A B : x ×m: (A + B) = x ×m: A + x ×m: B.

Lemma scalemxA x y A : x ×m: (y ×m: A) = (x × y) ×m: A.


Lemma scalemx_const a b : a *: const_mx b = const_mx (a × b).

Lemma matrix_sum_delta A :
  A = \sum_(i < m) \sum_(j < n) A i j *: delta_mx i j.

End RingModule.

Section StructuralLinear.

Lemma swizzle_mx_is_scalable m n p q f g k :
  scalable (@swizzle_mx R m n p q f g k).



Definition mxvec_is_linear m n := can2_scalable (@vec_mxK R m n) mxvecK.

End StructuralLinear.

Lemma trmx_delta m n i j : (delta_mx i j)^T = delta_mx j i :> 'M[R]_(n, m).

Lemma row_sum_delta n (u : 'rV_n) : u = \sum_(j < n) u 0 j *: delta_mx 0 j.

Lemma delta_mx_lshift m n1 n2 i j :
  delta_mx i (lshift n2 j) = row_mx (delta_mx i j) 0 :> 'M_(m, n1 + n2).

Lemma delta_mx_rshift m n1 n2 i j :
  delta_mx i (rshift n1 j) = row_mx 0 (delta_mx i j) :> 'M_(m, n1 + n2).

Lemma delta_mx_ushift m1 m2 n i j :
  delta_mx (lshift m2 i) j = col_mx (delta_mx i j) 0 :> 'M_(m1 + m2, n).

Lemma delta_mx_dshift m1 m2 n i j :
  delta_mx (rshift m1 i) j = col_mx 0 (delta_mx i j) :> 'M_(m1 + m2, n).

Lemma vec_mx_delta m n i j :
  vec_mx (delta_mx 0 (mxvec_index i j)) = delta_mx i j :> 'M_(m, n).

Lemma mxvec_delta m n i j :
  mxvec (delta_mx i j) = delta_mx 0 (mxvec_index i j) :> 'rV_(m × n).

Ltac split_mxE := apply/matrixPi j; do ![rewrite mxE | case: split ⇒ ?].

Lemma scale_row_mx m n1 n2 a (A1 : 'M_(m, n1)) (A2 : 'M_(m, n2)) :
  a *: row_mx A1 A2 = row_mx (a *: A1) (a *: A2).

Lemma scale_col_mx m1 m2 n a (A1 : 'M_(m1, n)) (A2 : 'M_(m2, n)) :
  a *: col_mx A1 A2 = col_mx (a *: A1) (a *: A2).

Lemma scale_block_mx m1 m2 n1 n2 a (Aul : 'M_(m1, n1)) (Aur : 'M_(m1, n2))
                                   (Adl : 'M_(m2, n1)) (Adr : 'M_(m2, n2)) :
  a *: block_mx Aul Aur Adl Adr
     = block_mx (a *: Aul) (a *: Aur) (a *: Adl) (a *: Adr).

Diagonal matrices
Scalar matrix : a diagonal matrix with a constant on the diagonal
Section ScalarMx.

Variable n : nat.

Fact scalar_mx_key : unit.
Definition scalar_mx x : 'M[R]_n :=
  \matrix[scalar_mx_key]_(i , j) (x *+ (i == j)).
Notation "x %:M" := (scalar_mx x) : ring_scope.

Lemma diag_const_mx a : diag_mx (const_mx a) = a%:M :> 'M_n.

Lemma tr_scalar_mx a : (a%:M)^T = a%:M.

Lemma trmx1 : (1%:M)^T = 1%:M.

Lemma scalar_mx_is_additive : additive scalar_mx.

Lemma scale_scalar_mx a1 a2 : a1 *: a2%:M = (a1 × a2)%:M :> 'M_n.

Lemma scalemx1 a : a *: 1%:M = a%:M.

Lemma scalar_mx_sum_delta a : a%:M = \sum_i a *: delta_mx i i.

Lemma mx1_sum_delta : 1%:M = \sum_i delta_mx i i.

Lemma row1 i : row i 1%:M = delta_mx 0 i.

Lemma col1 i : col i 1%:M = delta_mx i 0.

Definition is_scalar_mx (A : 'M[R]_n) :=
  if insub 0%N is Some i then A == (A i i)%:M else true.

Lemma is_scalar_mxP A : reflect ( a, A = a%:M) (is_scalar_mx A).

Lemma scalar_mx_is_scalar a : is_scalar_mx a%:M.

Lemma mx0_is_scalar : is_scalar_mx 0.

Lemma scalar_mx_is_diag a : is_diag_mx (a%:M).

Lemma is_scalar_mx_is_diag A : is_scalar_mx A is_diag_mx A.

Lemma scalar_mx_is_trig a : is_trig_mx (a%:M).

Lemma is_scalar_mx_is_trig A : is_scalar_mx A is_trig_mx A.

End ScalarMx.

Notation "x %:M" := (scalar_mx _ x) : ring_scope.

Lemma mx11_scalar (A : 'M_1) : A = (A 0 0)%:M.

Lemma scalar_mx_block n1 n2 a : a%:M = block_mx a%:M 0 0 a%:M :> 'M_(n1 + n2).

Matrix multiplication using bigops.
Fact mulmx_key : unit.
Definition mulmx {m n p} (A : 'M_(m, n)) (B : 'M_(n, p)) : 'M[R]_(m, p) :=
  \matrix[mulmx_key]_(i, k) \sum_j (A i j × B j k).


Lemma mulmxA m n p q (A : 'M_(m, n)) (B : 'M_(n, p)) (C : 'M_(p, q)) :
  A ×m (B ×m C) = A ×m B ×m C.

Lemma mul0mx m n p (A : 'M_(n, p)) : 0 ×m A = 0 :> 'M_(m, p).

Lemma mulmx0 m n p (A : 'M_(m, n)) : A ×m 0 = 0 :> 'M_(m, p).

Lemma mulmxN m n p (A : 'M_(m, n)) (B : 'M_(n, p)) : A ×m (- B) = - (A ×m B).

Lemma mulNmx m n p (A : 'M_(m, n)) (B : 'M_(n, p)) : - A ×m B = - (A ×m B).

Lemma mulmxDl m n p (A1 A2 : 'M_(m, n)) (B : 'M_(n, p)) :
  (A1 + A2) ×m B = A1 ×m B + A2 ×m B.

Lemma mulmxDr m n p (A : 'M_(m, n)) (B1 B2 : 'M_(n, p)) :
  A ×m (B1 + B2) = A ×m B1 + A ×m B2.

Lemma mulmxBl m n p (A1 A2 : 'M_(m, n)) (B : 'M_(n, p)) :
  (A1 - A2) ×m B = A1 ×m B - A2 ×m B.

Lemma mulmxBr m n p (A : 'M_(m, n)) (B1 B2 : 'M_(n, p)) :
  A ×m (B1 - B2) = A ×m B1 - A ×m B2.

Lemma mulmx_suml m n p (A : 'M_(n, p)) I r P (B_ : I 'M_(m, n)) :
   (\sum_(i <- r | P i) B_ i) ×m A = \sum_(i <- r | P i) B_ i ×m A.

Lemma mulmx_sumr m n p (A : 'M_(m, n)) I r P (B_ : I 'M_(n, p)) :
   A ×m (\sum_(i <- r | P i) B_ i) = \sum_(i <- r | P i) A ×m B_ i.

Lemma scalemxAl m n p a (A : 'M_(m, n)) (B : 'M_(n, p)) :
  a *: (A ×m B) = (a *: A) ×m B.
Right scaling associativity requires a commutative ring

Lemma rowE m n i (A : 'M_(m, n)) : row i A = delta_mx 0 i ×m A.

Lemma colE m n i (A : 'M_(m, n)) : col i A = A ×m delta_mx i 0.

Lemma mul_rVP m n A B :((@mulmx 1 m n)^~ A =1 mulmx^~ B) (A = B).

Lemma row_mul m n p (i : 'I_m) A (B : 'M_(n, p)) :
  row i (A ×m B) = row i A ×m B.

Lemma mulmx_sum_row m n (u : 'rV_m) (A : 'M_(m, n)) :
  u ×m A = \sum_i u 0 i *: row i A.

Lemma mxsub_mul m n m' n' p f g (A : 'M_(m, p)) (B : 'M_(p, n)) :
  mxsub f g (A ×m B) = rowsub f A ×m colsub g B :> 'M_(m', n').

Lemma mul_rowsub_mx m n m' p f (A : 'M_(m, p)) (B : 'M_(p, n)) :
  rowsub f A ×m B = rowsub f (A ×m B) :> 'M_(m', n).

Lemma mulmx_colsub m n n' p g (A : 'M_(m, p)) (B : 'M_(p, n)) :
  A ×m colsub g B = colsub g (A ×m B) :> 'M_(m, n').

Lemma mul_delta_mx_cond m n p (j1 j2 : 'I_n) (i1 : 'I_m) (k2 : 'I_p) :
  delta_mx i1 j1 ×m delta_mx j2 k2 = delta_mx i1 k2 *+ (j1 == j2).

Lemma mul_delta_mx m n p (j : 'I_n) (i : 'I_m) (k : 'I_p) :
  delta_mx i j ×m delta_mx j k = delta_mx i k.

Lemma mul_delta_mx_0 m n p (j1 j2 : 'I_n) (i1 : 'I_m) (k2 : 'I_p) :
  j1 != j2 delta_mx i1 j1 ×m delta_mx j2 k2 = 0.

Lemma mul_diag_mx m n d (A : 'M_(m, n)) :
  diag_mx d ×m A = \matrix_(i, j) (d 0 i × A i j).

Lemma mul_mx_diag m n (A : 'M_(m, n)) d :
  A ×m diag_mx d = \matrix_(i, j) (A i j × d 0 j).

Lemma mulmx_diag n (d e : 'rV_n) :
  diag_mx d ×m diag_mx e = diag_mx (\row_j (d 0 j × e 0 j)).

Lemma mul_scalar_mx m n a (A : 'M_(m, n)) : a%:M ×m A = a *: A.

Lemma scalar_mxM n a b : (a × b)%:M = a%:M ×m b%:M :> 'M_n.

Lemma mul1mx m n (A : 'M_(m, n)) : 1%:M ×m A = A.

Lemma mulmx1 m n (A : 'M_(m, n)) : A ×m 1%:M = A.

Lemma rowsubE m m' n f (A : 'M_(m, n)) :
   rowsub f A = rowsub f 1%:M ×m A :> 'M_(m', n).

mulmx and col_perm, row_perm, xcol, xrow

Lemma mul_col_perm m n p s (A : 'M_(m, n)) (B : 'M_(n, p)) :
  col_perm s A ×m B = A ×m row_perm s^-1 B.

Lemma mul_row_perm m n p s (A : 'M_(m, n)) (B : 'M_(n, p)) :
  A ×m row_perm s B = col_perm s^-1 A ×m B.

Lemma mul_xcol m n p j1 j2 (A : 'M_(m, n)) (B : 'M_(n, p)) :
  xcol j1 j2 A ×m B = A ×m xrow j1 j2 B.

Permutation matrix
Partial identity matrix (used in rank decomposition).
Block products; we cover all 1 x 2, 2 x 1, and 2 x 2 block products.
Lemma mul_mx_row m n p1 p2 (A : 'M_(m, n)) (Bl : 'M_(n, p1)) (Br : 'M_(n, p2)) :
  A ×m row_mx Bl Br = row_mx (A ×m Bl) (A ×m Br).

Lemma mul_col_mx m1 m2 n p (Au : 'M_(m1, n)) (Ad : 'M_(m2, n)) (B : 'M_(n, p)) :
  col_mx Au Ad ×m B = col_mx (Au ×m B) (Ad ×m B).

Lemma mul_row_col m n1 n2 p (Al : 'M_(m, n1)) (Ar : 'M_(m, n2))
                            (Bu : 'M_(n1, p)) (Bd : 'M_(n2, p)) :
  row_mx Al Ar ×m col_mx Bu Bd = Al ×m Bu + Ar ×m Bd.

Lemma mul_col_row m1 m2 n p1 p2 (Au : 'M_(m1, n)) (Ad : 'M_(m2, n))
                                (Bl : 'M_(n, p1)) (Br : 'M_(n, p2)) :
  col_mx Au Ad ×m row_mx Bl Br
     = block_mx (Au ×m Bl) (Au ×m Br) (Ad ×m Bl) (Ad ×m Br).

Lemma mul_row_block m n1 n2 p1 p2 (Al : 'M_(m, n1)) (Ar : 'M_(m, n2))
                                  (Bul : 'M_(n1, p1)) (Bur : 'M_(n1, p2))
                                  (Bdl : 'M_(n2, p1)) (Bdr : 'M_(n2, p2)) :
  row_mx Al Ar ×m block_mx Bul Bur Bdl Bdr
   = row_mx (Al ×m Bul + Ar ×m Bdl) (Al ×m Bur + Ar ×m Bdr).

Lemma mul_block_col m1 m2 n1 n2 p (Aul : 'M_(m1, n1)) (Aur : 'M_(m1, n2))
                                  (Adl : 'M_(m2, n1)) (Adr : 'M_(m2, n2))
                                  (Bu : 'M_(n1, p)) (Bd : 'M_(n2, p)) :
  block_mx Aul Aur Adl Adr ×m col_mx Bu Bd
   = col_mx (Aul ×m Bu + Aur ×m Bd) (Adl ×m Bu + Adr ×m Bd).

Lemma mulmx_block m1 m2 n1 n2 p1 p2 (Aul : 'M_(m1, n1)) (Aur : 'M_(m1, n2))
                                    (Adl : 'M_(m2, n1)) (Adr : 'M_(m2, n2))
                                    (Bul : 'M_(n1, p1)) (Bur : 'M_(n1, p2))
                                    (Bdl : 'M_(n2, p1)) (Bdr : 'M_(n2, p2)) :
  block_mx Aul Aur Adl Adr ×m block_mx Bul Bur Bdl Bdr
    = block_mx (Aul ×m Bul + Aur ×m Bdl) (Aul ×m Bur + Aur ×m Bdr)
               (Adl ×m Bul + Adr ×m Bdl) (Adl ×m Bur + Adr ×m Bdr).

Lemma mulmx_lsub m n p k (A : 'M_(m, n)) (B : 'M_(n, p + k)) :
  A ×m lsubmx B = lsubmx (A ×m B).

Lemma mulmx_rsub m n p k (A : 'M_(m, n)) (B : 'M_(n, p + k)) :
  A ×m rsubmx B = rsubmx (A ×m B).

Lemma mul_usub_mx m k n p (A : 'M_(m + k, n)) (B : 'M_(n, p)) :
  usubmx A ×m B = usubmx (A ×m B).

Lemma mul_dsub_mx m k n p (A : 'M_(m + k, n)) (B : 'M_(n, p)) :
  dsubmx A ×m B = dsubmx (A ×m B).

Correspondence between matrices and linear function on row vectors.
Section LinRowVector.

Variables m n : nat.

Fact lin1_mx_key : unit.
Definition lin1_mx (f : 'rV[R]_m 'rV[R]_n) :=
  \matrix[lin1_mx_key]_(i, j) f (delta_mx 0 i) 0 j.

Variable f : {linear 'rV[R]_m 'rV[R]_n}.

Lemma mul_rV_lin1 u : u ×m lin1_mx f = f u.

End LinRowVector.

Correspondence between matrices and linear function on matrices.
Section LinMatrix.

Variables m1 n1 m2 n2 : nat.

Definition lin_mx (f : 'M[R]_(m1, n1) 'M[R]_(m2, n2)) :=
  lin1_mx (mxvec \o f \o vec_mx).

Variable f : {linear 'M[R]_(m1, n1) 'M[R]_(m2, n2)}.

Lemma mul_rV_lin u : u ×m lin_mx f = mxvec (f (vec_mx u)).

Lemma mul_vec_lin A : mxvec A ×m lin_mx f = mxvec (f A).

Lemma mx_rV_lin u : vec_mx (u ×m lin_mx f) = f (vec_mx u).

Lemma mx_vec_lin A : vec_mx (mxvec A ×m lin_mx f) = f A.

End LinMatrix.

Section Mulmxr.

Variables m n p : nat.
Implicit Type A : 'M[R]_(m, n).
Implicit Type B : 'M[R]_(n, p).

Definition mulmxr B A := mulmx A B.
Arguments mulmxr B A /.

Definition lin_mulmxr B := lin_mx (mulmxr B).

Lemma mulmxr_is_linear B : linear (mulmxr B).

Lemma lin_mulmxr_is_linear : linear lin_mulmxr.

End Mulmxr.
Arguments mulmxr {_ _ _} B A /.

The trace.
Section Trace.

Variable n : nat.

Definition mxtrace (A : 'M[R]_n) := \sum_i A i i.

Lemma mxtrace_tr A : \tr A^T = \tr A.

Lemma mxtrace_is_scalar : scalar mxtrace.

Lemma mxtrace0 : \tr 0 = 0.
Lemma mxtraceD A B : \tr (A + B) = \tr A + \tr B.
Lemma mxtraceZ a A : \tr (a *: A) = a × \tr A.

Lemma mxtrace_diag D : \tr (diag_mx D) = \sum_j D 0 j.

Lemma mxtrace_scalar a : \tr a%:M = a *+ n.

Lemma mxtrace1 : \tr 1%:M = n%:R.

End Trace.

Lemma trace_mx11 (A : 'M_1) : \tr A = A 0 0.

Lemma mxtrace_block n1 n2 (Aul : 'M_n1) Aur Adl (Adr : 'M_n2) :
  \tr (block_mx Aul Aur Adl Adr) = \tr Aul + \tr Adr.

The matrix ring structure requires a strutural condition (dimension of the form n.+1) to satisfy the nontriviality condition we have imposed.
Section MatrixRing.

Variable n' : nat.

Lemma matrix_nonzero1 : 1%:M != 0 :> 'M_n.


Lemma mulmxE : mulmx = *%R.
Lemma idmxE : 1%:M = 1 :> 'M_n.

Lemma scalar_mx_is_multiplicative : multiplicative (@scalar_mx n).

End MatrixRing.

Section LiftPerm.

Block expresssion of a lifted permutation matrix, for the Cormen LUP.

Variable n : nat.

These could be in zmodp, but that would introduce a dependency on perm.

Definition lift0_perm s : 'S_n.+1 := lift_perm 0 0 s.

Lemma lift0_perm0 s : lift0_perm s 0 = 0.

Lemma lift0_perm_lift s k' :
  lift0_perm s (lift 0 k') = lift (0 : 'I_n.+1) (s k').

Lemma lift0_permK s : cancel (lift0_perm s) (lift0_perm s^-1).

Lemma lift0_perm_eq0 s i : (lift0_perm s i == 0) = (i == 0).

Block expresssion of a lifted permutation matrix

Definition lift0_mx A : 'M_(1 + n) := block_mx 1 0 0 A.

Lemma lift0_mx_perm s : lift0_mx (perm_mx s) = perm_mx (lift0_perm s).

Lemma lift0_mx_is_perm s : is_perm_mx (lift0_mx (perm_mx s)).

End LiftPerm.

Lemma exp_block_diag_mx m n (A: 'M_m.+1) (B : 'M_n.+1) k :
  (block_mx A 0 0 B) ^+ k = block_mx (A ^+ k) 0 0 (B ^+ k).

Determinants and adjugates are defined here, but most of their properties only hold for matrices over a commutative ring, so their theory is deferred to that section.
The determinant, in one line with the Leibniz Formula
Definition determinant n (A : 'M_n) : R :=
  \sum_(s : 'S_n) (-1) ^+ s × \prod_i A i (s i).

The cofactor of a matrix on the indexes i and j
Definition cofactor n A (i j : 'I_n) : R :=
  (-1) ^+ (i + j) × determinant (row' i (col' j A)).

The adjugate matrix : defined as the transpose of the matrix of cofactors
Fact adjugate_key : unit.
Definition adjugate n (A : 'M_n) := \matrix[adjugate_key]_(i, j) cofactor A j i.

End MatrixAlgebra.

Arguments delta_mx {R m n}.
Arguments scalar_mx {R n}.
Arguments perm_mx {R n}.
Arguments tperm_mx {R n}.
Arguments pid_mx {R m n}.
Arguments copid_mx {R n}.
Arguments lin_mulmxr {R m n p}.

Arguments is_scalar_mxP {R n A}.
Arguments mul_delta_mx {R m n p}.

#[global] Hint Extern 0 (is_true (is_diag_mx (scalar_mx _))) ⇒
  apply: scalar_mx_is_diag : core.
#[global] Hint Extern 0 (is_true (is_trig_mx (scalar_mx _))) ⇒
  apply: scalar_mx_is_trig : core.
#[global] Hint Extern 0 (is_true (is_diag_mx (diag_mx _))) ⇒
  apply: diag_mx_is_diag : core.
#[global] Hint Extern 0 (is_true (is_trig_mx (diag_mx _))) ⇒
  apply: diag_mx_is_trig : core.

Notation "a %:M" := (scalar_mx a) : ring_scope.
Notation "A *m B" := (mulmx A B) : ring_scope.
Arguments mulmxr {_ _ _ _} B A /.
Notation "\tr A" := (mxtrace A) : ring_scope.
Notation "'\det' A" := (determinant A) : ring_scope.
Notation "'\adj' A" := (adjugate A) : ring_scope.

Non-commutative transpose requires multiplication in the converse ring.
Lemma trmx_mul_rev (R : ringType) m n p (A : 'M[R]_(m, n)) (B : 'M[R]_(n, p)) :
  (A ×m B)^T = (B : 'M[R^c]_(n, p))^T ×m (A : 'M[R^c]_(m, n))^T.


Section FinZmodMatrix.
Variables (V : finZmodType) (m n : nat).


#[compress_coercions]
HB.instance Definition _ := [finGroupMixin of MV for +%R].

End FinZmodMatrix.

#[compress_coercions]
HB.instance Definition _ (R : finRingType) (m n : nat) :=
  FinRing.Zmodule.on 'M[R]_(m, n).

#[compress_coercions]
HB.instance Definition _ (R : finRingType) n :=
  [Finite of 'M[R]_n.+1 by <:].

Parametricity over the algebra structure.
Commutation property specialized to 'M[R]_n ************
GRing.comm is bound to (non trivial) rings, and matrices form a (non trivial) ring only when they are square and of manifestly positive size. However during proofs in endomorphism reduction, we take restrictions, which are matrices of size #|V| (with V a matrix space) and it becomes cumbersome to state commutation between restrictions, unless we relax the setting, and this relaxation corresponds to comm_mx A B := A *m B = B *m A. As witnessed by comm_mxE, when A and B have type 'M_n.+1, comm_mx A B is convertible to GRing.comm A B. The boolean version comm_mxb is designed to be used with seq.allrel

Context {R : ringType} {n : nat}.
Implicit Types (f g p : 'M[R]_n) (fs : seq 'M[R]_n) (d : 'rV[R]_n) (I : Type).

Definition comm_mx f g : Prop := f ×m g = g ×m f.
Definition comm_mxb f g : bool := f ×m g == g ×m f.

Lemma comm_mx_sym f g : comm_mx f g comm_mx g f.

Lemma comm_mx_refl f : comm_mx f f.

Lemma comm_mx0 f : comm_mx f 0.
Lemma comm0mx f : comm_mx 0 f.

Lemma comm_mx1 f : comm_mx f 1%:M.

Lemma comm1mx f : comm_mx 1%:M f.

Hint Resolve comm_mx0 comm0mx comm_mx1 comm1mx : core.

Lemma comm_mxN f g : comm_mx f g comm_mx f (- g).

Lemma comm_mxN1 f : comm_mx f (- 1%:M).

Lemma comm_mxD f g g' : comm_mx f g comm_mx f g' comm_mx f (g + g').

Lemma comm_mxB f g g' : comm_mx f g comm_mx f g' comm_mx f (g - g').

Lemma comm_mxM f g g' : comm_mx f g comm_mx f g' comm_mx f (g ×m g').

Lemma comm_mx_sum I (s : seq I) (P : pred I) (F : I 'M[R]_n) (f : 'M[R]_n) :
  ( i : I, P i comm_mx f (F i)) comm_mx f (\sum_(i <- s | P i) F i).

Lemma comm_mxP f g : reflect (comm_mx f g) (comm_mxb f g).

Notation all_comm_mx fs := (all2rel comm_mxb fs).

Lemma all_comm_mxP fs :
  reflect {in fs &, f g, f ×m g = g ×m f} (all_comm_mx fs).

Lemma all_comm_mx1 f : all_comm_mx [:: f].

Lemma all_comm_mx2P f g : reflect (f ×m g = g ×m f) (all_comm_mx [:: f; g]).

Lemma all_comm_mx_cons f fs :
  all_comm_mx (f :: fs) = all (comm_mxb f) fs && all_comm_mx fs.

End CommMx.
Notation all_comm_mx := (allrel comm_mxb).

Lemma comm_mxE (R : ringType) (n : nat) : @comm_mx R n.+1 = @GRing.comm _.

Section ComMatrix.
Lemmas for matrices with coefficients in a commutative ring
Variable R : comRingType.

Section AssocLeft.

Variables m n p : nat.
Implicit Type A : 'M[R]_(m, n).
Implicit Type B : 'M[R]_(n, p).

Lemma trmx_mul A B : (A ×m B)^T = B^T ×m A^T.

Lemma scalemxAr a A B : a *: (A ×m B) = A ×m (a *: B).

Lemma mulmx_is_scalable A : scalable (@mulmx _ m n p A).

Definition lin_mulmx A : 'M[R]_(n × p, m × p) := lin_mx (mulmx A).

Lemma lin_mulmx_is_linear : linear lin_mulmx.

End AssocLeft.

Section LinMulRow.

Variables m n : nat.

Definition lin_mul_row u : 'M[R]_(m × n, n) := lin1_mx (mulmx u \o vec_mx).

Lemma lin_mul_row_is_linear : linear lin_mul_row.

Lemma mul_vec_lin_row A u : mxvec A ×m lin_mul_row u = u ×m A.

End LinMulRow.

Lemma mxvec_dotmul m n (A : 'M[R]_(m, n)) u v :
  mxvec (u^T ×m v) ×m (mxvec A)^T = u ×m A ×m v^T.

Section MatrixAlgType.

Variable n' : nat.


End MatrixAlgType.

Lemma diag_mxC n (d e : 'rV[R]_n) :
  diag_mx d ×m diag_mx e = diag_mx e ×m diag_mx d.

Lemma diag_mx_comm n (d e : 'rV[R]_n) : comm_mx (diag_mx d) (diag_mx e).

Lemma scalar_mxC m n a (A : 'M[R]_(m, n)) : A ×m a%:M = a%:M ×m A.

Lemma mul_mx_scalar m n a (A : 'M[R]_(m, n)) : A ×m a%:M = a *: A.

Lemma comm_mx_scalar n a (A : 'M[R]_n) : comm_mx A a%:M.

Lemma comm_scalar_mx n a (A : 'M[R]_n) : comm_mx a%:M A.

Lemma mxtrace_mulC m n (A : 'M[R]_(m, n)) B :
   \tr (A ×m B) = \tr (B ×m A).

The theory of determinants

Lemma determinant_multilinear n (A B C : 'M[R]_n) i0 b c :
    row i0 A = b *: row i0 B + c *: row i0 C
    row' i0 B = row' i0 A
    row' i0 C = row' i0 A
  \det A = b × \det B + c × \det C.

Lemma determinant_alternate n (A : 'M[R]_n) i1 i2 :
  i1 != i2 A i1 =1 A i2 \det A = 0.

Lemma det_tr n (A : 'M[R]_n) : \det A^T = \det A.

Lemma det_perm n (s : 'S_n) : \det (perm_mx s) = (-1) ^+ s :> R.

Lemma det1 n : \det (1%:M : 'M[R]_n) = 1.

Lemma det_mx00 (A : 'M[R]_0) : \det A = 1.

Lemma detZ n a (A : 'M[R]_n) : \det (a *: A) = a ^+ n × \det A.

Lemma det0 n' : \det (0 : 'M[R]_n'.+1) = 0.

Lemma det_scalar n a : \det (a%:M : 'M[R]_n) = a ^+ n.

Lemma det_scalar1 a : \det (a%:M : 'M[R]_1) = a.

Lemma det_mx11 (M : 'M[R]_1) : \det M = M 0 0.

Lemma det_mulmx n (A B : 'M[R]_n) : \det (A ×m B) = \det A × \det B.

Lemma detM n' (A B : 'M[R]_n'.+1) : \det (A × B) = \det A × \det B.

Laplace expansion lemma
Lemma expand_cofactor n (A : 'M[R]_n) i j :
  cofactor A i j =
    \sum_(s : 'S_n | s i == j) (-1) ^+ s × \prod_(k | i != k) A k (s k).

Lemma expand_det_row n (A : 'M[R]_n) i0 :
  \det A = \sum_j A i0 j × cofactor A i0 j.

Lemma cofactor_tr n (A : 'M[R]_n) i j : cofactor A^T i j = cofactor A j i.

Lemma cofactorZ n a (A : 'M[R]_n) i j :
  cofactor (a *: A) i j = a ^+ n.-1 × cofactor A i j.

Lemma expand_det_col n (A : 'M[R]_n) j0 :
  \det A = \sum_i (A i j0 × cofactor A i j0).

Lemma trmx_adj n (A : 'M[R]_n) : (\adj A)^T = \adj A^T.

Lemma adjZ n a (A : 'M[R]_n) : \adj (a *: A) = a^+n.-1 *: \adj A.

Cramer Rule : adjugate on the left
Cramer rule : adjugate on the right
Lemma mul_adj_mx n (A : 'M[R]_n) : \adj A ×m A = (\det A)%:M.

Lemma adj1 n : \adj (1%:M) = 1%:M :> 'M[R]_n.

Left inverses are right inverses.
Lemma mulmx1C n (A B : 'M[R]_n) : A ×m B = 1%:M B ×m A = 1%:M.

Only tall matrices have inverses.
Lemma mulmx1_min m n (A : 'M[R]_(m, n)) B : A ×m B = 1%:M m n.

Lemma det_ublock n1 n2 Aul (Aur : 'M[R]_(n1, n2)) Adr :
  \det (block_mx Aul Aur 0 Adr) = \det Aul × \det Adr.

Lemma det_lblock n1 n2 Aul (Adl : 'M[R]_(n2, n1)) Adr :
  \det (block_mx Aul 0 Adl Adr) = \det Aul × \det Adr.

Lemma det_trig n (A : 'M[R]_n) : is_trig_mx A \det A = \prod_(i < n) A i i.

Lemma det_diag n (d : 'rV[R]_n) : \det (diag_mx d) = \prod_i d 0 i.

End ComMatrix.

Arguments lin_mul_row {R m n} u.
Arguments lin_mulmx {R m n p} A.
Arguments comm_mx_scalar {R n}.
Arguments comm_scalar_mx {R n}.
Arguments diag_mx_comm {R n}.


#[global] Hint Resolve comm_mx_scalar comm_scalar_mx : core.

Matrix unit ring and inverse matrices ***************

Section MatrixInv.

Variables R : comUnitRingType.

Section Defs.

Variable n : nat.
Implicit Type A : 'M[R]_n.

Definition unitmx : pred 'M[R]_n := fun A\det A \is a GRing.unit.
Definition invmx A := if A \in unitmx then (\det A)^-1 *: \adj A else A.

Lemma unitmxE A : (A \in unitmx) = (\det A \is a GRing.unit).

Lemma unitmx1 : 1%:M \in unitmx.

Lemma unitmx_perm s : perm_mx s \in unitmx.

Lemma unitmx_tr A : (A^T \in unitmx) = (A \in unitmx).

Lemma unitmxZ a A : a \is a GRing.unit (a *: A \in unitmx) = (A \in unitmx).

Lemma invmx1 : invmx 1%:M = 1%:M.

Lemma invmxZ a A : a *: A \in unitmx invmx (a *: A) = a^-1 *: invmx A.

Lemma invmx_scalar a : invmx (a%:M) = a^-1%:M.

Lemma mulVmx : {in unitmx, left_inverse 1%:M invmx mulmx}.

Lemma mulmxV : {in unitmx, right_inverse 1%:M invmx mulmx}.

Lemma mulKmx m : {in unitmx, @left_loop _ 'M_(n, m) invmx mulmx}.

Lemma mulKVmx m : {in unitmx, @rev_left_loop _ 'M_(n, m) invmx mulmx}.

Lemma mulmxK m : {in unitmx, @right_loop 'M_(m, n) _ invmx mulmx}.

Lemma mulmxKV m : {in unitmx, @rev_right_loop 'M_(m, n) _ invmx mulmx}.

Lemma det_inv A : \det (invmx A) = (\det A)^-1.

Lemma unitmx_inv A : (invmx A \in unitmx) = (A \in unitmx).

Lemma unitmx_mul A B : (A ×m B \in unitmx) = (A \in unitmx) && (B \in unitmx).

Lemma trmx_inv (A : 'M_n) : (invmx A)^T = invmx (A^T).

Lemma invmxK : involutive invmx.

Lemma mulmx1_unit A B : A ×m B = 1%:M A \in unitmx B \in unitmx.

Lemma intro_unitmx A B : B ×m A = 1%:M A ×m B = 1%:M unitmx A.

Lemma invmx_out : {in [predC unitmx], invmx =1 id}.

End Defs.

Variable n' : nat.


Lemmas requiring that the coefficients are in a unit ring
Finite inversible matrices and the general linear group.
Section FinUnitMatrix.

Variables (n : nat) (R : finComUnitRingType).

Definition GLtype of phant R := {unit 'M[R]_n.-1.+1}.

Coercion GLval ph (u : GLtype ph) : 'M[R]_n.-1.+1 :=
  let: FinRing.Unit A _ := u in A.

End FinUnitMatrix.

Bind Scope group_scope with GLtype.
Arguments GLval {n%N R ph} u%g.

Notation "{ ''GL_' n [ R ] }" := (GLtype n (Phant R))
  (at level 0, n at level 2, format "{ ''GL_' n [ R ] }") : type_scope.
Notation "{ ''GL_' n ( p ) }" := {'GL_n['F_p]}
  (at level 0, n at level 2, p at level 10,
    format "{ ''GL_' n ( p ) }") : type_scope.


Section GL_unit.

Variables (n : nat) (R : finComUnitRingType).


Definition GLgroup of phant R := [set: {'GL_n[R]}].
Canonical GLgroup_group ph := Eval hnf in [group of GLgroup ph].

Implicit Types u v : {'GL_n[R]}.

Lemma GL_1E : GLval 1 = 1.
Lemma GL_VE u : GLval u^-1 = (GLval u)^-1.
Lemma GL_VxE u : GLval u^-1 = invmx u.
Lemma GL_ME u v : GLval (u × v) = GLval u × GLval v.
Lemma GL_MxE u v : GLval (u × v) = u ×m v.
Lemma GL_unit u : GLval u \is a GRing.unit.
Lemma GL_unitmx u : val u \in unitmx.

Lemma GL_det u : \det u != 0.

End GL_unit.

Notation "''GL_' n [ R ]" := (GLgroup n (Phant R))
  (at level 8, n at level 2, format "''GL_' n [ R ]") : group_scope.
Notation "''GL_' n ( p )" := 'GL_n['F_p]
  (at level 8, n at level 2, p at level 10,
   format "''GL_' n ( p )") : group_scope.
Notation "''GL_' n [ R ]" := (GLgroup_group n (Phant R)) : Group_scope.
Notation "''GL_' n ( p )" := (GLgroup_group n (Phant 'F_p)) : Group_scope.

Matrices over a domain ******************************

Section MatrixDomain.

Variable R : idomainType.

Lemma scalemx_eq0 m n a (A : 'M[R]_(m, n)) :
  (a *: A == 0) = (a == 0) || (A == 0).

Lemma scalemx_inj m n a :
  a != 0 injective ( *:%R a : 'M[R]_(m, n) 'M[R]_(m, n)).

Lemma det0P n (A : 'M[R]_n) :
  reflect (exists2 v : 'rV[R]_n, v != 0 & v ×m A = 0) (\det A == 0).

End MatrixDomain.

Arguments det0P {R n A}.

Parametricity at the field level (mx_is_scalar, unit and inverse are only mapped at this level).
Section MapFieldMatrix.

Variables (aF : fieldType) (rF : comUnitRingType) (f : {rmorphism aF rF}).

Lemma map_mx_inj {m n} : injective (map_mx f : 'M_(m, n) 'M_(m, n)).

Lemma map_mx_is_scalar n (A : 'M_n) : is_scalar_mx A^f = is_scalar_mx A.

Lemma map_unitmx n (A : 'M_n) : (A^f \in unitmx) = (A \in unitmx).

Lemma map_mx_unit n' (A : 'M_n'.+1) :
  (A^f \is a GRing.unit) = (A \is a GRing.unit).

Lemma map_invmx n (A : 'M_n) : (invmx A)^f = invmx A^f.

Lemma map_mx_inv n' (A : 'M_n'.+1) : A^-1^f = A^f^-1.

Lemma map_mx_eq0 m n (A : 'M_(m, n)) : (A^f == 0) = (A == 0).

End MapFieldMatrix.

Arguments map_mx_inj {aF rF f m n} [A1 A2] eqA12f : rename.

LUP decomposion *****************************

Section CormenLUP.

Variable F : fieldType.

Decomposition of the matrix A to P A = L U with
  • P a permutation matrix
  • L a unipotent lower triangular matrix
  • U an upper triangular matrix

Fixpoint cormen_lup {n} :=
  match n return let M := 'M[F]_n.+1 in M M × M × M with
  | 0 ⇒ fun A(1, 1, A)
  | _.+1fun A
    let k := odflt 0 [pick k | A k 0 != 0] in
    let A1 : 'M_(1 + _) := xrow 0 k A in
    let P1 : 'M_(1 + _) := tperm_mx 0 k in
    let Schur := ((A k 0)^-1 *: dlsubmx A1) ×m ursubmx A1 in
    let: (P2, L2, U2) := cormen_lup (drsubmx A1 - Schur) in
    let P := block_mx 1 0 0 P2 ×m P1 in
    let L := block_mx 1 0 ((A k 0)^-1 *: (P2 ×m dlsubmx A1)) L2 in
    let U := block_mx (ulsubmx A1) (ursubmx A1) 0 U2 in
    (P, L, U)
  end.

Lemma cormen_lup_perm n (A : 'M_n.+1) : is_perm_mx (cormen_lup A).1.1.

Lemma cormen_lup_correct n (A : 'M_n.+1) :
  let: (P, L, U) := cormen_lup A in P × A = L × U.

Lemma cormen_lup_detL n (A : 'M_n.+1) : \det (cormen_lup A).1.2 = 1.

Lemma cormen_lup_lower n A (i j : 'I_n.+1) :
  i j (cormen_lup A).1.2 i j = (i == j)%:R.

Lemma cormen_lup_upper n A (i j : 'I_n.+1) :
  j < i (cormen_lup A).2 i j = 0 :> F.

End CormenLUP.

Section mxOver.
Section mxOverType.
Context {m n : nat} {T : Type}.
Implicit Types (S : {pred T}).

Definition mxOver_pred (S : {pred T}) :=
  fun M : 'M[T]_(m, n)[ i, [ j, M i j \in S]].
Arguments mxOver_pred _ _ /.
Definition mxOver (S : {pred T}) := [qualify a M | mxOver_pred S M].

Lemma mxOverP {S : {pred T}} {M : 'M[T]__} :
  reflect ( i j, M i j \in S) (M \is a mxOver S).

Lemma mxOverS (S1 S2 : {pred T}) :
  {subset S1 S2} {subset mxOver S1 mxOver S2}.

Lemma mxOver_const c S : c \in S const_mx c \is a mxOver S.

Lemma mxOver_constE c S : (m > 0)%N (n > 0)%N
  (const_mx c \is a mxOver S) = (c \in S).

End mxOverType.

Lemma thinmxOver {n : nat} {T : Type} (M : 'M[T]_(n, 0)) S : M \is a mxOver S.

Lemma flatmxOver {n : nat} {T : Type} (M : 'M[T]_(0, n)) S : M \is a mxOver S.

Section mxOverZmodule.
Context {M : zmodType} {m n : nat}.
Implicit Types (S : {pred M}).

Lemma mxOver0 S : 0 \in S 0 \is a @mxOver m n _ S.

Section mxOverAdd.
Variable addS : addrClosed M.
Fact mxOver_add_subproof : addr_closed (@mxOver m n _ addS).
End mxOverAdd.

Section mxOverOpp.
Variable oppS : opprClosed M.
Fact mxOver_opp_subproof : oppr_closed (@mxOver m n _ oppS).
End mxOverOpp.


End mxOverZmodule.

Section mxOverRing.
Context {R : ringType} {m n : nat}.

Lemma mxOver_scalar S c : 0 \in S c \in S c%:M \is a @mxOver n n R S.

Lemma mxOver_scalarE S c : (n > 0)%N
  (c%:M \is a @mxOver n n R S) = ((n > 1) ==> (0 \in S)) && (c \in S).

Lemma mxOverZ (S : mulrClosed R) :
  {in S & mxOver S, a : R, v : 'M[R]_(m, n),
     a *: v \is a mxOver S}.

Lemma mxOver_diag (S : {pred R}) k (D : 'rV[R]_k) :
   0 \in S D \is a mxOver S diag_mx D \is a mxOver S.

Lemma mxOver_diagE (S : {pred R}) k (D : 'rV[R]_k) : k > 0
  (diag_mx D \is a mxOver S) = ((k > 1) ==> (0 \in S)) && (D \is a mxOver S).

Lemma mxOverM (S : semiringClosed R) p q r : {in mxOver S & mxOver S,
   u : 'M[R]_(p, q), v : 'M[R]_(q, r), u ×m v \is a mxOver S}.

End mxOverRing.

Section mxRingOver.
Context {R : ringType} {n : nat}.

Section semiring.
Variable S : semiringClosed R.
Fact mxOver_mul_subproof : mulr_closed (@mxOver n.+1 n.+1 _ S).
End semiring.


End mxRingOver.
End mxOver.

Section BlockMatrix.
Import tagnat.
Context {T : Type} {p q : nat} {p_ : 'I_p nat} {q_ : 'I_q nat}.
Notation sp := (\sum_i p_ i)%N.
Notation sq := (\sum_i q_ i)%N.
Implicit Type (s : 'I_sp) (t : 'I_sq).

Definition mxblock (B_ : i j, 'M[T]_(p_ i, q_ j)) :=
  \matrix_(j, k) B_ (sig1 j) (sig1 k) (sig2 j) (sig2 k).

Definition mxrow m (B_ : j, 'M[T]_(m, q_ j)) :=
  \matrix_(j, k) B_ (sig1 k) j (sig2 k).

Definition mxcol n (B_ : i, 'M[T]_(p_ i, n)) :=
  \matrix_(j, k) B_ (sig1 j) (sig2 j) k.

Definition submxblock (A : 'M[T]_(sp, sq)) i j := mxsub (Rank i) (Rank j) A.
Definition submxrow m (A : 'M[T]_(m, sq)) j := colsub (Rank j) A.
Definition submxcol n (A : 'M[T]_(sp, n)) i := rowsub (Rank i) A.

Lemma mxblockEh B_ : \mxblock_(i, j) B_ i j = \mxrow_j \mxcol_i B_ i j.

Lemma mxblockEv B_ : \mxblock_(i, j) B_ i j = \mxcol_i \mxrow_j B_ i j.

Lemma submxblockEh A i j : submxblock A i j = submxcol (submxrow A j) i.

Lemma submxblockEv A i j : submxblock A i j = submxrow (submxcol A i) j.

Lemma mxblockK B_ i j : submxblock (\mxblock_(i, j) B_ i j) i j = B_ i j.

Lemma mxrowK m B_ j : @submxrow m (\mxrow_j B_ j) j = B_ j.

Lemma mxcolK n B_ i : @submxcol n (\mxcol_i B_ i) i = B_ i.

Lemma submxrow_matrix B_ j :
  submxrow (\mxblock_(i, j) B_ i j) j = \mxcol_i B_ i j.

Lemma submxcol_matrix B_ i :
  submxcol (\mxblock_(i, j) B_ i j) i = \mxrow_j B_ i j.

Lemma submxblockK A : \mxblock_(i, j) (submxblock A i j) = A.

Lemma submxrowK m (A : 'M[T]_(m, sq)) : \mxrow_j (submxrow A j) = A.

Lemma submxcolK n (A : 'M[T]_(sp, n)) : \mxcol_i (submxcol A i) = A.

Lemma mxblockP A B :
  ( i j, submxblock A i j = submxblock B i j) A = B.

Lemma mxrowP m (A B : 'M_(m, sq)) :
  ( j, submxrow A j = submxrow B j) A = B.

Lemma mxcolP n (A B : 'M_(sp, n)) :
  ( i, submxcol A i = submxcol B i) A = B.

Lemma eq_mxblockP A_ B_ :
  ( i j, A_ i j = B_ i j)
  (\mxblock_(i, j) A_ i j = \mxblock_(i, j) B_ i j).

Lemma eq_mxblock A_ B_ :
  ( i j, A_ i j = B_ i j)
  (\mxblock_(i, j) A_ i j = \mxblock_(i, j) B_ i j).

Lemma eq_mxrowP m (A_ B_ : j, 'M[T]_(m, q_ j)) :
  ( j, A_ j = B_ j) (\mxrow_j A_ j = \mxrow_j B_ j).

Lemma eq_mxrow m (A_ B_ : j, 'M[T]_(m, q_ j)) :
  ( j, A_ j = B_ j) (\mxrow_j A_ j = \mxrow_j B_ j).

Lemma eq_mxcolP n (A_ B_ : i, 'M[T]_(p_ i, n)) :
  ( i, A_ i = B_ i) (\mxcol_i A_ i = \mxcol_i B_ i).

Lemma eq_mxcol n (A_ B_ : i, 'M[T]_(p_ i, n)) :
  ( i, A_ i = B_ i) (\mxcol_i A_ i = \mxcol_i B_ i).

Lemma row_mxrow m (B_ : j, 'M[T]_(m, q_ j)) i :
  row i (\mxrow_j B_ j) = \mxrow_j (row i (B_ j)).

Lemma col_mxrow m (B_ : j, 'M[T]_(m, q_ j)) j :
  col j (\mxrow_j B_ j) = col (sig2 j) (B_ (sig1 j)).

Lemma row_mxcol n (B_ : i, 'M[T]_(p_ i, n)) i :
  row i (\mxcol_i B_ i) = row (sig2 i) (B_ (sig1 i)).

Lemma col_mxcol n (B_ : i, 'M[T]_(p_ i, n)) j :
  col j (\mxcol_i B_ i) = \mxcol_i (col j (B_ i)).

Lemma row_mxblock B_ i :
  row i (\mxblock_(i, j) B_ i j) = \mxrow_j row (sig2 i) (B_ (sig1 i) j).

Lemma col_mxblock B_ j :
  col j (\mxblock_(i, j) B_ i j) = \mxcol_i col (sig2 j) (B_ i (sig1 j)).

End BlockMatrix.

Notation "\mxblock_ ( i < m , j < n ) E" :=
  (mxblock (fun (i : 'I_m) (j : 'I_ n) ⇒ E)) (only parsing) : ring_scope.
Notation "\mxblock_ ( i , j < n ) E" :=
  (\mxblock_(i < n, j < n) E) (only parsing) : ring_scope.
Notation "\mxblock_ ( i , j ) E" := (\mxblock_(i < _, j < _) E) : ring_scope.
Notation "\mxrow_ ( j < m ) E" := (mxrow (fun (j : 'I_m) ⇒ E))
  (only parsing) : ring_scope.
Notation "\mxrow_ j E" := (\mxrow_(j < _) E) : ring_scope.
Notation "\mxcol_ ( i < m ) E" := (mxcol (fun (i : 'I_m) ⇒ E))
  (only parsing) : ring_scope.
Notation "\mxcol_ i E" := (\mxcol_(i < _) E) : ring_scope.

Lemma tr_mxblock {T : Type} {p q : nat} {p_ : 'I_p nat} {q_ : 'I_q nat}
  (B_ : i j, 'M[T]_(p_ i, q_ j)) :
  (\mxblock_(i, j) B_ i j)^T = \mxblock_(i, j) (B_ j i)^T.

Section SquareBlockMatrix.

Context {T : Type} {p : nat} {p_ : 'I_p nat}.
Notation sp := (\sum_i p_ i)%N.
Implicit Type (s : 'I_sp).

Lemma tr_mxrow n (B_ : j, 'M[T]_(n, p_ j)) :
  (\mxrow_j B_ j)^T = \mxcol_i (B_ i)^T.

Lemma tr_mxcol n (B_ : i, 'M[T]_(p_ i, n)) :
  (\mxcol_i B_ i)^T = \mxrow_i (B_ i)^T.

Lemma tr_submxblock (A : 'M[T]_sp) i j :
  (submxblock A i j)^T = (submxblock A^T j i).

Lemma tr_submxrow n (A : 'M[T]_(n, sp)) j :
  (submxrow A j)^T = (submxcol A^T j).

Lemma tr_submxcol n (A : 'M[T]_(sp, n)) i :
  (submxcol A i)^T = (submxrow A^T i).

End SquareBlockMatrix.

Section BlockRowRecL.
Import tagnat.
Context {T : Type} {m : nat} {p_ : 'I_m.+1 nat}.
Notation sp := (\sum_i p_ i)%N.

Lemma mxsize_recl : (p_ ord0 + \sum_i p_ (lift ord0 i) = (\sum_i p_ i))%N.

Lemma mxrow_recl n (B_ : j, 'M[T]_(n, p_ j)) :
  \mxrow_j B_ j = castmx (erefl, mxsize_recl)
    (row_mx (B_ 0) (\mxrow_j B_ (lift ord0 j))).

End BlockRowRecL.

Lemma mxcol_recu {T : Type} {p : nat} {p_ : 'I_p.+1 nat} m
    (B_ : j, 'M[T]_(p_ j, m)) :
  \mxcol_j B_ j = castmx (mxsize_recl, erefl)
    (col_mx (B_ 0) (\mxcol_j B_ (lift ord0 j))).

Section BlockMatrixRec.
Context {T : Type}.

Lemma mxblock_recu {p q : nat} {p_ : 'I_p.+1 nat} {q_ : 'I_q nat}
    (B_ : i j, 'M[T]_(p_ i, q_ j)) :
  \mxblock_(i, j) B_ i j = castmx (mxsize_recl, erefl) (col_mx
     (\mxrow_j B_ ord0 j)
     (\mxblock_(i, j) B_ (l0 i) j)).

Lemma mxblock_recl {p q : nat} {p_ : 'I_p nat} {q_ : 'I_q.+1 nat}
    (B_ : i j, 'M[T]_(p_ i, q_ j)) :
  \mxblock_(i, j) B_ i j = castmx (erefl, mxsize_recl)
  (row_mx (\mxcol_i B_ i ord0) (\mxblock_(i, j) B_ i (l0 j))).

Lemma mxblock_recul {p q : nat} {p_ : 'I_p.+1 nat} {q_ : 'I_q.+1 nat}
    (B_ : i j, 'M[T]_(p_ i, q_ j)) :
  \mxblock_(i, j) B_ i j = castmx e (block_mx
     (B_ 0 0) (\mxrow_j B_ ord0 (l0 j))
     (\mxcol_i B_ (l0 i) ord0) (\mxblock_(i, j) B_ (l0 i) (l0 j))).

Lemma mxrowEblock {q : nat} {q_ : 'I_q nat} m
    (R_ : j, 'M[T]_(m, q_ j)) :
  (\mxrow_j R_ j) =
  castmx (big_ord1 _ (fun m), erefl) (\mxblock_(i < 1, j < q) R_ j).

Lemma mxcolEblock {p : nat} {p_ : 'I_p nat} n
    (C_ : i, 'M[T]_(p_ i, n)) :
  (\mxcol_i C_ i) =
  castmx (erefl, big_ord1 _ (fun n)) (\mxblock_(i < p, j < 1) C_ i).

Lemma mxEmxrow m n (A : 'M[T]_(m, n)) :
  A = castmx (erefl, big_ord1 _ (fun n)) (\mxrow__ A).

Lemma mxEmxcol m n (A : 'M[T]_(m, n)) :
  A = castmx (big_ord1 _ (fun m), erefl) (\mxcol__ A).

Lemma mxEmxblock m n (A : 'M[T]_(m, n)) :
  A = castmx (big_ord1 _ (fun m), big_ord1 _ (fun n))
             (\mxblock_(i < 1, j < 1) A).

End BlockMatrixRec.

Section BlockRowZmod.
Context {V : zmodType} {q : nat} {q_ : 'I_q nat}.
Notation sq := (\sum_i q_ i)%N.
Implicit Type (s : 'I_sq).

Lemma mxrowD m (R_ R'_ : j, 'M[V]_(m, q_ j)) :
  \mxrow_j (R_ j + R'_ j) = \mxrow_j (R_ j) + \mxrow_j (R'_ j).

Lemma mxrowN m (R_ : j, 'M[V]_(m, q_ j)) :
  \mxrow_j (- R_ j) = - \mxrow_j (R_ j).

Lemma mxrowB m (R_ R'_ : j, 'M[V]_(m, q_ j)) :
  \mxrow_j (R_ j - R'_ j) = \mxrow_j (R_ j) - \mxrow_j (R'_ j).

Lemma mxrow0 m : \mxrow_j (0 : 'M[V]_(m, q_ j)) = 0.

Lemma mxrow_const m a : \mxrow_j (const_mx a : 'M[V]_(m, q_ j)) = const_mx a.

Lemma mxrow_sum (J : finType) m
    (R_ : i j, 'M[V]_(m, q_ j)) (P : {pred J}) :
  \mxrow_j (\sum_(i | P i) R_ i j) = \sum_(i | P i) \mxrow_j (R_ i j).

Lemma submxrowD m (B B' : 'M[V]_(m, sq)) j :
 submxrow (B + B') j = submxrow B j + submxrow B' j.

Lemma submxrowN m (B : 'M[V]_(m, sq)) j :
 submxrow (- B) j = - submxrow B j.

Lemma submxrowB m (B B' : 'M[V]_(m, sq)) j :
 submxrow (B - B') j = submxrow B j - submxrow B' j.

Lemma submxrow0 m j : submxrow (0 : 'M[V]_(m, sq)) j = 0.

Lemma submxrow_sum (J : finType) m
   (R_ : i, 'M[V]_(m, sq)) (P : {pred J}) j:
  submxrow (\sum_(i | P i) R_ i) j = \sum_(i | P i) submxrow (R_ i) j.

End BlockRowZmod.

Section BlockRowRing.
Context {R : ringType} {n : nat} {q_ : 'I_n nat}.
Notation sq := (\sum_i q_ i)%N.
Implicit Type (s : 'I_sq).

Lemma mul_mxrow m n' (A : 'M[R]_(m, n')) (R_ : j, 'M[R]_(n', q_ j)) :
  A ×m \mxrow_j R_ j= \mxrow_j (A ×m R_ j).

Lemma mul_submxrow m n' (A : 'M[R]_(m, n')) (B : 'M[R]_(n', sq)) j :
  A ×m submxrow B j= submxrow (A ×m B) j.

End BlockRowRing.

Section BlockColZmod.
Context {V : zmodType} {n : nat} {p_ : 'I_n nat}.
Notation sp := (\sum_i p_ i)%N.
Implicit Type (s : 'I_sp).

Lemma mxcolD m (C_ C'_ : i, 'M[V]_(p_ i, m)) :
  \mxcol_i (C_ i + C'_ i) = \mxcol_i (C_ i) + \mxcol_i (C'_ i).

Lemma mxcolN m (C_ : i, 'M[V]_(p_ i, m)) :
  \mxcol_i (- C_ i) = - \mxcol_i (C_ i).

Lemma mxcolB m (C_ C'_ : i, 'M[V]_(p_ i, m)) :
  \mxcol_i (C_ i - C'_ i) = \mxcol_i (C_ i) - \mxcol_i (C'_ i).

Lemma mxcol0 m : \mxcol_i (0 : 'M[V]_(p_ i, m)) = 0.

Lemma mxcol_const m a : \mxcol_j (const_mx a : 'M[V]_(p_ j, m)) = const_mx a.

Lemma mxcol_sum (I : finType) m (C_ : j i, 'M[V]_(p_ i, m)) (P : {pred I}):
  \mxcol_i (\sum_(j | P j) C_ j i) = \sum_(j | P j) \mxcol_i (C_ j i).

Lemma submxcolD m (B B' : 'M[V]_(sp, m)) i :
 submxcol (B + B') i = submxcol B i + submxcol B' i.

Lemma submxcolN m (B : 'M[V]_(sp, m)) i :
 submxcol (- B) i = - submxcol B i.

Lemma submxcolB m (B B' : 'M[V]_(sp, m)) i :
 submxcol (B - B') i = submxcol B i - submxcol B' i.

Lemma submxcol0 m i : submxcol (0 : 'M[V]_(sp, m)) i = 0.

Lemma submxcol_sum (I : finType) m
   (C_ : j, 'M[V]_(sp, m)) (P : {pred I}) i :
  submxcol (\sum_(j | P j) C_ j) i = \sum_(j | P j) submxcol (C_ j) i.

End BlockColZmod.

Section BlockColRing.
Context {R : ringType} {n : nat} {p_ : 'I_n nat}.
Notation sp := (\sum_i p_ i)%N.
Implicit Type (s : 'I_sp).

Lemma mxcol_mul n' m (C_ : i, 'M[R]_(p_ i, n')) (A : 'M[R]_(n', m)) :
  \mxcol_i C_ i ×m A = \mxcol_i (C_ i ×m A).

Lemma submxcol_mul n' m (B : 'M[R]_(sp, n')) (A : 'M[R]_(n', m)) i :
  submxcol B i ×m A = submxcol (B ×m A) i.

End BlockColRing.

Section BlockMatrixZmod.
Context {V : zmodType} {m n : nat}.
Context {p_ : 'I_m nat} {q_ : 'I_n nat}.
Notation sp := (\sum_i p_ i)%N.
Notation sq := (\sum_i q_ i)%N.

Lemma mxblockD (B_ B'_ : i j, 'M[V]_(p_ i, q_ j)) :
  \mxblock_(i, j) (B_ i j + B'_ i j) =
  \mxblock_(i, j) (B_ i j) + \mxblock_(i, j) (B'_ i j).

Lemma mxblockN (B_ : i j, 'M[V]_(p_ i, q_ j)) :
  \mxblock_(i, j) (- B_ i j) = - \mxblock_(i, j) (B_ i j).

Lemma mxblockB (B_ B'_ : i j, 'M[V]_(p_ i, q_ j)) :
  \mxblock_(i, j) (B_ i j - B'_ i j) =
  \mxblock_(i, j) (B_ i j) - \mxblock_(i, j) (B'_ i j).

Lemma mxblock0 : \mxblock_(i, j) (0 : 'M[V]_(p_ i, q_ j)) = 0.

Lemma mxblock_const a :
  \mxblock_(i, j) (const_mx a : 'M[V]_(p_ i, q_ j)) = const_mx a.

Lemma mxblock_sum (I : finType)
    (B_ : k i j, 'M[V]_(p_ i, q_ j)) (P : {pred I}):
  \mxblock_(i, j) (\sum_(k | P k) B_ k i j) =
  \sum_(k | P k) \mxblock_(i, j) (B_ k i j).

Lemma submxblockD (B B' : 'M[V]_(sp, sq)) i j :
 submxblock (B + B') i j = submxblock B i j + submxblock B' i j.

Lemma submxblockN (B : 'M[V]_(sp, sq)) i j :
 submxblock (- B) i j = - submxblock B i j.

Lemma submxblockB (B B' : 'M[V]_(sp, sq)) i j :
 submxblock (B - B') i j = submxblock B i j - submxblock B' i j.

Lemma submxblock0 i j : submxblock (0 : 'M[V]_(sp, sq)) i j = 0.

Lemma submxblock_sum (I : finType)
   (B_ : k, 'M[V]_(sp, sq)) (P : {pred I}) i j :
  submxblock (\sum_(k | P k) B_ k) i j = \sum_(k | P k) submxblock (B_ k) i j.

End BlockMatrixZmod.

Section BlockMatrixRing.
Context {R : ringType} {p q : nat} {p_ : 'I_p nat} {q_ : 'I_q nat}.
Notation sp := (\sum_i p_ i)%N.
Notation sq := (\sum_i q_ i)%N.

Lemma mul_mxrow_mxcol m n
    (R_ : j, 'M[R]_(m, p_ j)) (C_ : i, 'M[R]_(p_ i, n)) :
  \mxrow_j R_ j ×m \mxcol_i C_ i = \sum_i (R_ i ×m C_ i).

Lemma mul_mxcol_mxrow m
    (C_ : i, 'M[R]_(p_ i, m)) (R_ : j, 'M[R]_(m, q_ j)) :
  \mxcol_i C_ i×m \mxrow_j R_ j = \mxblock_(i, j) (C_ i ×m R_ j).

Lemma mul_mxrow_mxblock m
    (R_ : i, 'M[R]_(m, p_ i)) (B_ : i j, 'M[R]_(p_ i, q_ j)) :
  \mxrow_i R_ i ×m \mxblock_(i, j) B_ i j = \mxrow_j (\sum_i (R_ i ×m B_ i j)).

Lemma mul_mxblock_mxrow m
    (B_ : i j, 'M[R]_(q_ i, p_ j)) (C_ : i, 'M[R]_(p_ i, m)) :
  \mxblock_(i, j) B_ i j ×m \mxcol_j C_ j = \mxcol_i (\sum_j (B_ i j ×m C_ j)).

End BlockMatrixRing.

Lemma mul_mxblock {R : ringType} {p q r : nat}
    {p_ : 'I_p nat} {q_ : 'I_q nat} {r_ : 'I_r nat}
    (A_ : i j, 'M[R]_(p_ i, q_ j)) (B_ : j k, 'M_(q_ j, r_ k)) :
  \mxblock_(i, j) A_ i j ×m \mxblock_(j, k) B_ j k =
  \mxblock_(i, k) \sum_j (A_ i j ×m B_ j k).

Section SquareBlockMatrixZmod.
Import Order.TTheory tagnat.
Context {V : zmodType} {p : nat} {p_ : 'I_p nat}.
Notation sp := (\sum_i p_ i)%N.
Implicit Type (s : 'I_sp).

Lemma is_trig_mxblockP (B_ : i j, 'M[V]_(p_ i, p_ j)) :
  reflect [/\ (i j : 'I_p), (i < j)%N B_ i j = 0 &
               i, is_trig_mx (B_ i i)]
          (is_trig_mx (\mxblock_(i, j) B_ i j)).

Lemma is_trig_mxblock (B_ : i j, 'M[V]_(p_ i, p_ j)) :
  is_trig_mx (\mxblock_(i, j) B_ i j) =
  ([ i : 'I_p, j : 'I_p, (i < j)%N ==> (B_ i j == 0)] &&
   [ i, is_trig_mx (B_ i i)]).

Lemma is_diag_mxblockP (B_ : i j, 'M[V]_(p_ i, p_ j)) :
  reflect [/\ (i j : 'I_p), i != j B_ i j = 0 &
               i, is_diag_mx (B_ i i)]
          (is_diag_mx (\mxblock_(i, j) B_ i j)).

Lemma is_diag_mxblock (B_ : i j, 'M[V]_(p_ i, p_ j)) :
  is_diag_mx (\mxblock_(i, j) B_ i j) =
  ([ i : 'I_p, j : 'I_p, (i != j) ==> (B_ i j == 0)] &&
   [ i, is_diag_mx (B_ i i)]).

Definition mxdiag (B_ : i, 'M[V]_(p_ i)) : 'M[V]_(\sum_i p_ i) :=
  \mxblock_(j, k) if j == k then conform_mx 0 (B_ j) else 0.

Lemma submxblock_diag (B_ : i, 'M[V]_(p_ i)) i :
  submxblock (\mxdiag_i B_ i) i i = B_ i.

Lemma eq_mxdiagP (B_ B'_ : i, 'M[V]_(p_ i)) :
  ( i, B_ i = B'_ i) (\mxdiag_i B_ i = \mxdiag_i B'_ i).

Lemma eq_mxdiag (B_ B'_ : i, 'M[V]_(p_ i)) :
  ( i, B_ i = B'_ i) (\mxdiag_i B_ i = \mxdiag_i B'_ i).

Lemma mxdiagD (B_ B'_ : i, 'M[V]_(p_ i)) :
  \mxdiag_i (B_ i + B'_ i) = \mxdiag_i (B_ i) + \mxdiag_i (B'_ i).

Lemma mxdiagN (B_ : i, 'M[V]_(p_ i)) :
  \mxdiag_i (- B_ i) = - \mxdiag_i (B_ i).

Lemma mxdiagB (B_ B'_ : i, 'M[V]_(p_ i)) :
  \mxdiag_i (B_ i - B'_ i) = \mxdiag_i (B_ i) - \mxdiag_i (B'_ i).

Lemma mxdiag0 : \mxdiag_i (0 : 'M[V]_(p_ i)) = 0.

Lemma mxdiag_sum (I : finType) (B_ : k i, 'M[V]_(p_ i)) (P : {pred I}) :
  \mxdiag_i (\sum_(k | P k) B_ k i) = \sum_(k | P k) \mxdiag_i (B_ k i).

Lemma tr_mxdiag (B_ : i, 'M[V]_(p_ i)) :
  (\mxdiag_i B_ i)^T = \mxdiag_i (B_ i)^T.

Lemma row_mxdiag (B_ : i, 'M[V]_(p_ i)) k :
  let B'_ i := if sig1 k == i then conform_mx 0 (B_ i) else 0 in
  row k (\mxdiag_ i B_ i) = row (sig2 k) (\mxrow_i B'_ i).

Lemma col_mxdiag (B_ : i, 'M[V]_(p_ i)) k :
  let B'_ i := if sig1 k == i then conform_mx 0 (B_ i) else 0 in
  col k (\mxdiag_ i B_ i) = col (sig2 k) (\mxcol_i B'_ i).

End SquareBlockMatrixZmod.

Notation "\mxdiag_ ( i < n ) E" := (mxdiag (fun i : 'I_nE))
  (only parsing) : ring_scope.
Notation "\mxdiag_ i E" := (\mxdiag_(i < _) E) : ring_scope.

Lemma mxdiag_recl {V : zmodType} {m : nat} {p_ : 'I_m.+1 nat}
    (B_ : i, 'M[V]_(p_ i)) :
  \mxdiag_i B_ i = castmx (mxsize_recl, mxsize_recl)
    (block_mx (B_ 0) 0 0 (\mxdiag_i B_ (lift ord0 i))).

Section SquareBlockMatrixRing.
Import tagnat.
Context {R : ringType} {p : nat} {p_ : 'I_p nat}.
Notation sp := (\sum_i p_ i)%N.
Implicit Type (s : 'I_sp).

Lemma mxtrace_mxblock (B_ : i j, 'M[R]_(p_ i, p_ j)) :
  \tr (\mxblock_(i, j) B_ i j) = \sum_i \tr (B_ i i).

Lemma mxdiagZ a : \mxdiag_i (a%:M : 'M[R]_(p_ i)) = a%:M.

Lemma diag_mxrow (B_ : j, 'rV[R]_(p_ j)) :
  diag_mx (\mxrow_j B_ j) = \mxdiag_j (diag_mx (B_ j)).

Lemma mxtrace_mxdiag (B_ : i, 'M[R]_(p_ i)) :
  \tr (\mxdiag_i B_ i) = \sum_i \tr (B_ i).

Lemma mul_mxdiag_mxcol m
    (D_ : i, 'M[R]_(p_ i)) (C_ : i, 'M[R]_(p_ i, m)):
  \mxdiag_i D_ i ×m \mxcol_i C_ i = \mxcol_i (D_ i ×m C_ i).

End SquareBlockMatrixRing.

Lemma mul_mxrow_mxdiag {R : ringType} {p : nat} {p_ : 'I_p nat} m
    (R_ : i, 'M[R]_(m, p_ i)) (D_ : i, 'M[R]_(p_ i)) :
  \mxrow_i R_ i ×m \mxdiag_i D_ i = \mxrow_i (R_ i ×m D_ i).

Lemma mul_mxblock_mxdiag {R : ringType} {p q : nat}
  {p_ : 'I_p nat} {q_ : 'I_q nat}
    (B_ : i j, 'M[R]_(p_ i, q_ j)) (D_ : j, 'M[R]_(q_ j)) :
  \mxblock_(i, j) B_ i j ×m \mxdiag_j D_ j = \mxblock_(i, j) (B_ i j ×m D_ j).

Lemma mul_mxdiag_mxblock {R : ringType} {p q : nat}
  {p_ : 'I_p nat} {q_ : 'I_q nat}
    (D_ : j, 'M[R]_(p_ j)) (B_ : i j, 'M[R]_(p_ i, q_ j)):
  \mxdiag_j D_ j ×m \mxblock_(i, j) B_ i j = \mxblock_(i, j) (D_ i ×m B_ i j).