Add LoadPath "./interval".

Require Import Fourier Rdefinitions Raxioms RIneq Rbasic_fun Rfunctions 
  R_sqrt.
Require Import  ssreflect ssrbool ssrnat eqtype ssrfun seq fintype finfun 
   finset matrix bigops ssralg perm.
Require Import Rstruct min_max topology interval_R 
  real_matrix eigen intMatrix.

Set Implicit Arguments.
Unset Strict Implicit.
Import Prenex Implicits.

Open Scope R_scope.
Delimit Scope R_scope with Re.

Open Scope ring_scope.
Delimit Scope ring_scope with Ri.

Import GRing.Theory.

Section MatrixIntSys.

Variable m n: nat.
Variable A: matrix closedInt m n.
Variable b: matrix closedInt m 1.

(* need to coerce interval matrices and vectors to sets (preds) 
  of matrices and vectors of reals or do it manually?*)

Theorem Beeck1:
  forall x, sigma_sol A b x <-> exists t , setI (inSetm (mmuls_i A x)) (inSetm b) t.
Proof.
move=> x; split=> [Hsol|Hinter].
  elim Hsol=> vA [Ha [vB [Hb Hx]]].
  exists vB;  split; last by done.
  by rewrite -Hx; apply mmuls_i_correct.
elim Hinter => vb [Hvb1 Hvb2].
rewrite /sigma_sol.
elim ( mvmuls_i_compl Hvb1) => vA [Ha1 Ha2] .
exists vA; split; first by done.
exists vb; split; first by done.
by rewrite Ha2.
Qed.


Theorem Beeck2:
 forall x, sigma_sol A b x <-> inSetm (madd_i (mmuls_i A x) (mopp_i b)) 0 .
Proof.
move=> x; rewrite Beeck1; split => Hs.
  elim Hs => t [HtA Htb].
  replace (GRing.zero (matrix_zmodType real_zmodType m (S O)))
    with ( t + (- t)); [|rewrite -matrixP=> i j; rewrite !mxE; toR; ring].
  apply madd_i_correct; first by done.
  by apply mopp_i_correct.
  elim (madd_i_compl Hs) => [av  [Hav1 [Hav2 Havplus]]].
  exists av.1; split; first by done.
  have Hop: av.2 = - av.1.
  move: Havplus;  rewrite -!matrixP => Havplus i j; have hav:= Havplus i j.
  rewrite !mxE; rewrite !mxE in hav.
  replace (@GRing.zero real_zmodType) with (av.1 i j + (-av.1 i j))%Re in hav; 
    last by (toR; ring).
  by symmetry; apply Rplus_eq_reg_l with (av.1 i j).
rewrite Hop in Hav2.
replace b with (mopp_i (mopp_i b)).
  replace (av.1) with (- (- av.1)).
    by apply mopp_i_correct.
  by rewrite  -matrixP => i j; rewrite !mxE; rewrite ->Ropp_involutive.
rewrite /mopp_i -matrixP=> i j; rewrite !mxE /opp_i eq_intervalP //=.
by split; rewrite Ropp_involutive.
Qed.
   

Corollary OetteliPrager:
 forall x, sigma_sol A b x <-> 
   Mabs (mmid A *m x  - mmid b) <=m:
        mrad A *m Mabs x + mrad b.
Proof.
move=> x.
rewrite Beeck1.
have Ha: forall u v, 
  (exists t :  'cV_m, setI (inSetm u) (inSetm v) t) <->
  Mabs (mmid u - mmid v) <=m: mrad u + mrad v.
move=> u v; split => H1.
elim H1=> t [Htu Htv].
move=> i j; rewrite !mxE // .
have Heq: (mid (u i j) - mid (v i j))%Ri = 
  ((mid (u i j) - t i j) + (t i j - mid (v i j)))%Re by toR; ring.
rewrite Heq.
apply Rle_trans with  (1:=Rabs_triang _ _ ).
apply Rplus_le_compat.
by rewrite Rabs_minus_sym -in_mid_rad.
by rewrite -in_mid_rad.
have H3:  forall i, Rabs (mid (u i ord0) - mid (v i ord0)) 
  <= rad (u i ord0) + rad (v i ord0) ->
     {r : R | r \in (u i ord0) /\ r \in (v i ord0)}.
by move=>*; apply: intersect_int.
have H4:  forall i : 'I_m,
     Rabs (mid (u i ord0) - mid (v i ord0)) <= rad (u i ord0) + rad (v i ord0).
move=> i; have H11:= H1 i ord0.
by rewrite !mxE in H11.
have H5: forall i, {r : R | r \in u i ord0 /\ r \in v i ord0}.
by move=> i; apply H3.
exists (\col_i sval (H5 i)).
by split; move=> i j; rewrite mxE ; case: (H5 i) => [r [Hr1 Hr2]]; 
  rewrite //= (ord1_0 j).
set a:= (mmuls_i A x).
have Hmid: mmid a = mmid A *m x.
rewrite /a -matrixP=> i k; rewrite !mxE mid_big_add.
apply eq_bigr=> j _.
by rewrite /mmid mxE -mid_mul.
have Hrad: mrad a = mrad A *m Mabs x.
rewrite /a -matrixP=> i k; rewrite !mxE rad_big_add.
apply eq_bigr=> j _.
by rewrite /mrad !mxE  -rad_mul.
rewrite -Hmid -Hrad.
apply: Ha.
Qed.

End MatrixIntSys.


Section basicRegCrit.

Variable n': nat.
Notation n := n'.+1. 
Variable A: matrix closedInt n n.


Lemma reg_all_sol0:
  regular A <-> forall (x: 'cV_n) a, inSetm A a ->  a *m x = 0 -> x = 0.
Proof.
split => H; [move=> x a Ha Hax|].
have H0:= H _ Ha.
have H1: (invmx a) *m  (a *m x) =  (invmx a) *m 0.
by congr (_ *m _).
rewrite mulmxA in H1.
have H2: (invmx a *m a) = 1%:M.
rewrite -mulVmx; toR.
rewrite H2 mul1mx in H1.
by rewrite H1 mulmx0.
move=> a Ha.
have H1 := fun x => H x a Ha.
have Hdet:= det0P_cv_n (a: matrix real_field _ _).
by apply/Hdet.
Qed.


Lemma reg_int_sol0:
  regular A <-> (forall x: 'cV_n, inSetm (mmuls_i A x) 0 -> x = 0).
 rewrite reg_all_sol0.
split => H x.
  move=> Hin; elim : (mvmuls_i_compl Hin) => a [Ha1 Ha2].
  by apply: (H _ a); [| rewrite Ha2].
by move=> a Ha Hax; apply: H; rewrite -Hax; apply: mmuls_i_correct.
Qed.


Lemma reg_sig_sol0:
regular A <-> 
(forall x, sigma_sol A (\col_i (R0:closedInt)) x <-> x = 0).
Proof.
rewrite reg_all_sol0.
split=> H x.
  split=> Hss.
    elim: Hss => a [Ha [b [Hb Hab]]].
    apply: (H _ a); first by done.
    rewrite Hab -matrixP => i k; have Hbi := Hb i k.
    move: Hbi; rewrite in_int !mxE //= => Hbi.
    by apply Rle_antisym; elim Hbi.
  exists (minf A); split.
    by move=> i j; rewrite mxE //=; apply: inf_in_int.
  exists (\col_i R0); split.
    by move=> i k; rewrite !mxE in_int //=; split; right.
  by rewrite Hss mulmx0.
move=> a Ha Hax; rewrite -H; exists a; split; first by done.
exists 0; split; last by done.
by move=> i k; rewrite !mxE in_int; split; right.
Qed.

Lemma reg_ineq_mid_rad:
  regular A <-> 
  (forall x: 'cV_n, Mabs (mmid A *m x) <=m:  mrad A *m Mabs x -> x =  0).
Proof.
rewrite reg_sig_sol0; split=> H x; have Hx:= H x.
  move: Hx; rewrite OetteliPrager => Hx.
  move=> Hv; rewrite -Hx => i k; rewrite !mxE /Rminus.
  apply Rle_trans with (1:= Rabs_triang _ _).
  apply Rplus_le_compat.
    by have Hvi:= (Hv i k); rewrite !mxE in Hvi.
  right; rewrite Rabs_Ropp /mid /rad /= Rabs_right; [field|fourier].
split.
  rewrite OetteliPrager => Hop; apply Hx=> i k; rewrite !mxE.
  have Hopi:= (Hop i k); rewrite !mxE in Hopi.
  apply Rle_trans with  (Rabs (\sum_j mmid A i j * x j k - mid R0)).
    right; f_equal; rewrite /mid /=; field.
  apply Rle_trans with (1:= Hopi).
    right; f_equal; rewrite /rad /=; toR; field.
move=> Hx0.
exists (minf A); split.
  by move=> i j; rewrite mxE //=; apply: inf_in_int.
exists (\col_i R0); split.
  by move=> i k; rewrite !mxE in_int //=; split; right.
by rewrite Hx0 mulmx0.
Qed.


(* classicaly equivalent to previous lemma *)
Require Import Classical.

Lemma sing_ineq_mid_rad:
singular A <-> 
  (exists x: 'cV_n, Mabs (mmid A *m x) <=m: (mrad A *m Mabs x) /\ x <> 0).
Proof.
rewrite iff_not -reg_nsing_eq  reg_ineq_mid_rad.
split => H.
apply all_not_not_ex=> vA.
apply or_not_and.
apply imply_to_or => Ha Heq; elim Heq.
by apply H.
apply: not_ex_not_all => H1; elim H1=> x Hx.
apply H; exists x.
by apply imply_to_and.
Qed.

Lemma sing_ineq_mid_rad_norm:
singular A <-> 
  (exists x: 'cV_n, Mabs (mmid A *m x) <=m: (mrad A *m Mabs x) /\ x^T *m x = 1).
Proof.
rewrite sing_ineq_mid_rad; split; move=> [x [Hineq Hx]].
exists ((1/ sqrt ((x^T *m x) ord0 ord0)) *m: x).
split.
rewrite Mabs_scalemx -!scalemxAr Mabs_scalemx -!mul_scalar_mx.
apply Mmul_mle_compat_l; first by done.
move=> i j; rewrite !mxE; case C: (i == j); rewrite //=; 
  [rewrite mulr1n; apply Rabs_pos | by rewrite mulr0n; right].
rewrite trmx_scale -scalemxAr -scalemxAl scalemxA.
move/eqP: Hx => Hx.
move/vec_prod_pos: Hx => Hx.
have Hx1:= Hx ord0 ord0.
rewrite mxE in Hx1.
have Hsqrt:= sqrt_lt_R0 _ Hx1.

have Heq:  ((1 / sqrt ((x^T *m x) ord0 ord0) * (1 / sqrt ((x^T *m x) ord0 ord0))))
  = (1/ (x^T *m x) ord0 ord0).
toR; ring_simplify.
have Csqrt: (sqrt ((x^T *m x) ord0 ord0) != 0).
by apply/eqP; apply Rgt_not_eq.
have Cx: (x^T *m x) ord0 ord0 != 0.
by apply/eqP; apply Rgt_not_eq.
rewrite /Rinvx Csqrt Cx -Rinv_pow; last by apply Rgt_not_eq. 
by rewrite /pow Rmult_1_r sqrt_sqrt; last by left.
rewrite Heq.
apply/matrixP => i j; rewrite (ord1_0 i) (ord1_0 j) mxE -mulrA mulVr.
by rewrite  mxE /= mulr1n mulr1 .
by rewrite /GRing.unit //= /unit_R; apply/eqP; apply Rgt_not_eq.

exists x; split; first by done.
apply/eqP; apply/vec_prod_pos; rewrite Hx => i j; 
  rewrite !mxE (ord1_0 i) (ord1_0 j) /= mulr1n; apply Rlt_0_1.
Qed.


Lemma oettli_prager_sing:
  {x: 'cV_n |  Mabs ( mmid A *m x) <=m:  mrad A *m Mabs x /\ 
  x <> 0} ->   {vA | inSetm A vA /\ \det vA == 0}.
Proof.
move=>[x [Hre Hneq]].
  set y := (\col_i if ( (mrad A *m Mabs x) i ord0 != 0) then 
                                (mmid A *m x) i ord0 / ((mrad A) *m (Mabs x)) i ord0
                            else 1).
  pose z := \col_(i<n) (if (Rleb 0 (x i ord0)) then R1 else (-R1)).
  pose vA := \matrix_(i, j) ((mmid A) i j - y i ord0 * z j ord0 *  (mrad A) i j).
  exists vA; split.
    move=> i j; rewrite in_mid_rad.
    rewrite /vA mxE.
    replace ( (mmid A i j - y i ord0 * z j ord0 * mrad A i j)%Ri - mid (A i j))%Re
  (*  replace ((mid (A i j) + -( y i ord0 * z j ord0 * rad (A i j)) - mid (A i j)))%Re *)
  with
      ((- y i ord0 * z j ord0 * rad (A i j))%Re). (*; [|toR; ring].*)
    rewrite Rabs_mult -{2}(Rmult_1_l (rad (A i j))) (Rabs_right (rad _));
    [| apply Rle_ge; apply rad_pos].
    apply Rmult_le_compat_r; first by apply rad_pos .
    rewrite /z /y mxE; 
      rewrite (mxE  ( fun i0 _=> (if Rleb 0 (x i0 ord0) then 1 else -1)) j ).
    case Cle : (Rleb 0 (x j ord0)); rewrite Rabs_mult.
       rewrite (Rabs_right 1); [|fourier].
       rewrite  Rmult_1_r.
       case Cle2: ( (mrad A *m Mabs x) i ord0 != 0); 
         [|rewrite Rabs_Ropp Rabs_right; toR; fourier].
       rewrite Rabs_Ropp /Rdiv Rabs_mult.
         apply Rmult_le_reg_l with (Rabs ((mrad A *m Mabs x) i ord0 )).
         by apply Rabs_pos_lt; move/eqP: Cle2.
      rewrite Rmult_comm Rmult_assoc -Rabs_mult.
      have Htor:= @GRing.mulVr real_unitRing ( (mrad A *m Mabs x) i ord0).
      rewrite /GRing.mul /GRing.inv /cancel /= in Htor; toR; rewrite Htor;
      last by apply Cle2.
      rewrite Rabs_R1 !Rmult_1_r.
      have Hrei:= Hre i ord0.
      rewrite mxE in Hrei; apply Rle_trans with (1:=Hrei).
      rewrite Rabs_right; first by right.
      rewrite mxE.
      apply big_prop; first by right.
        by move=> a b Hz Ht; apply Rle_ge; apply: Rplus_le_le_0_compat; 
        apply Rge_le.
      move=> k _; rewrite mxE mxE; apply Rle_ge; apply Rmult_pos;
      [apply rad_pos| apply Rabs_pos].
    rewrite !Rabs_Ropp Rabs_R1 Rmult_1_r.
    case Cle2: ( (mrad A *m Mabs x) i ord0  != 0);
         [|rewrite Rabs_R1; right; ring].
  rewrite /Rdiv Rabs_mult.
         apply Rmult_le_reg_l with (Rabs ((mrad A *m Mabs x) i ord0)).
         by apply Rabs_pos_lt; move/eqP: Cle2.
      rewrite Rmult_comm Rmult_assoc -Rabs_mult.
      have Htor:= @GRing.mulVr real_unitRing ( (mrad A *m Mabs x) i ord0).
      rewrite /GRing.mul /GRing.inv /cancel /= in Htor; toR; rewrite Htor;
      last by apply Cle2.
      rewrite Rabs_R1 !Rmult_1_r.
      have Hrei:= Hre i ord0.
      rewrite mxE in Hrei; apply Rle_trans with (1:=Hrei).
      rewrite Rabs_right; first by right.
      rewrite mxE.
      apply big_prop; first by right.
        by move=> a b Hz Ht; apply Rle_ge; apply: Rplus_le_le_0_compat; 
        apply Rge_le.
      move=> k _; rewrite mxE mxE; apply Rle_ge; apply Rmult_pos;
      [apply rad_pos| apply Rabs_pos].
rewrite /mmid /mrad !mxE; toR; ring.
have Hdet:= det0P_cv_n (vA: matrix real_field _ _).
have Heq:  (\det vA == 0) =  (~~ (\det vA != 0)).
case ((\det vA == 0)); by rewrite negbK.

(*by split=> H; apply/negPn.
apply/eqP.*)
rewrite Heq.
apply/Hdet => H1.
apply Hneq.
apply H1.


(*apply/negP.
move/eqP => H1.*)
(*elim Hneq.
 rewrite Hdet => H1; apply H1.*)
apply/matrixP => i k; rewrite !mxE /vA .
transitivity (\sum_j (mmid A i j - y i ord0 * z j ord0 * mrad A i j)* x j ord0 ).
by apply eq_bigr => j _;  rewrite  mxE (ord1_0 k).
transitivity 
  (\sum_j ((mmid A i j * x j ord0 ) - (y i ord0 * (z j ord0 * mrad A i j * x j ord0)) )).
apply eq_bigr => j _; by toR; ring.
rewrite GRing.sumr_sub -big_distrr /=.
rewrite /y mxE;  case Cle: ((mrad A *m Mabs x) i ord0 != 0).
rewrite !mxE.
have Heq1:  (\sum_j mrad A i j * Mabs x j ord0) = 
  (\sum_(i0 < n) z i0 ord0 * mrad A i i0 * x i0 ord0) .
apply eq_bigr => j _; rewrite (GRing.mulrC (z j ord0)) -GRing.mulrA.
apply Rmult_eq_compat_l.
rewrite /Mabs  /z !mxE /Rabs.
case Cle2: (Rleb 0 (x j ord0)); move/RleP: Cle2 => Cle2.
by case: (Rcase_abs (x j ord0)) => Cle3; [elim (Rle_not_lt _ _ Cle2)| toR; ring].
by case:  (Rcase_abs (x j ord0)) => Cle3; [toR; ring| elim Cle2; apply Rge_le].
rewrite -Heq1 -mulrA mulrC  mulVr; 
  last by (rewrite mxE in Cle; apply Cle).
rewrite mul1r; toR; ring.
move/eqP: Cle => Cle.
rewrite mxE in Cle.
have Heq1:  (\sum_j mrad A i j * Mabs x j ord0) = 
  (\sum_(i0 < n) z i0 ord0 * mrad A i i0 * x i0 ord0) .
apply eq_bigr => j _; rewrite (GRing.mulrC (z j ord0)) -GRing.mulrA.
apply Rmult_eq_compat_l.
rewrite /Mabs  /z !mxE /Rabs.
case Cle2: (Rleb 0 (x j ord0)); move/RleP: Cle2 => Cle2.
by case:  (Rcase_abs (x j ord0)) => Cle3; [elim (Rle_not_lt _ _ Cle2)| toR;ring].
by case:  (Rcase_abs (x j ord0)) => Cle3; [toR; ring| elim Cle2; apply Rge_le].
rewrite -Heq1 Cle mulr0 oppr0 addr0.
have Hrei := Hre i ord0; rewrite !mxE Cle in Hrei.
have H: Rabs (\sum_j mmid A i j * x j ord0) = 0.
by apply Rle_antisym; [|apply Rabs_pos].
move: H; rewrite /Rabs;  case (Rcase_abs ) => Ca H .
rewrite -{2}oppr0 in H.
rewrite <-(Ropp_involutive _).
rewrite <-(Ropp_involutive ).
apply Ropp_eq_compat; rewrite H; by toR.
done.
Qed.


End basicRegCrit.


Section midInv_RegCrit.

Variable n': nat.
Notation n:= n'.+1.

Variable A: matrix closedInt n n.

(* uses proof by contradiction *)
Theorem thm31_midinva_specr:
forall es X, 
(forall M lam, lam \in (eigenv M) <-> (lam \in (es M))) ->
specr es (Mabs (1 - X * mmid A) + Mabs X * mrad A )  < 1 -> regular A.
Proof.
move=> es M Hlam Hsr.
apply nsing_reg; rewrite sing_ineq_mid_rad => Hsing.
elim Hsing => x [Hle Hdif0].

have Hinq: Mabs x <=m:  
  (Mabs (1 - M * mmid A) + Mabs M * mrad A )  *m (Mabs x).
have H1:  x =   (1 - M * mmid A) *m x +  (M * mmid A) *m x.
rewrite mulmx_subl mul1mx -matrixP => i k ; rewrite !mxE; toR; ring.
rewrite {1}H1.
apply Mle_trans with (1:= Mabs_triang _ _).
apply Mle_trans with 
( Mabs (1 - M * mmid A) *m (Mabs  x) +  
   (Mabs M ) *m (Mabs (mmid A *m x))) .
apply Madd_mle_compat; [apply Mabs_mul|].
rewrite -mulmxA; apply Mabs_mul.
rewrite mulmx_addl; apply Madd_mle_compat.
by move=> i j; right.
rewrite -mulmxA; apply Mmul_mle_compat_l; first by done.
move=> i j; rewrite !mxE; apply Rabs_pos.
have Heq:  1 *m: (Mabs x ) = Mabs x.
by apply/matrixP=> i j; rewrite mxE; toR; ring.
move: Hinq; rewrite -{1}Heq => Hinq.
pose a := (Mabs (1 - M * mmid A) + Mabs M * mrad A) .
rewrite -/a in Hsr Hinq.
move : Hsr; apply Rle_not_lt; rewrite corPF; last by done.
by exists (Mabs x); split; [move=> i j; rewrite !mxE; apply Rabs_pos|].
rewrite /a; move => i j; rewrite !mxE.
apply Rplus_le_le_0_compat; [apply Rabs_pos|].
apply: big_prop; [by right| apply Rplus_le_le_0_compat|].
move=> k _; rewrite !mxE; apply Rmult_pos; [apply Rabs_pos| apply rad_pos].
Qed.



Corollary cor32_midinv_specr:
forall es, \det (mmid A) != 0 -> 
(forall M lam, lam \in (eigenv M) <-> (lam \in (es M))) ->
  specr es (Mabs (mmid A)^-1 * mrad A )  < 1 -> regular A.
Proof.
move=> es Hma Hlam Hsr.
apply (@thm31_midinva_specr es ( (mmid A)^-1)).
done.
rewrite GRing.mulVr.
rewrite  GRing.addrN.
have Heq : @Mabs n n 0 = 0.
rewrite /Mabs ; apply/matrixP => i j; rewrite !mxE; apply Rabs_R0.
by rewrite Heq GRing.add0r.
by apply Hma.
Qed.


Theorem thm33:
(exists Rm , exists j,  
   \col_i (1 + Mabs ( 1 - mmid A * Rm))%Ri i j <=m:  
   \col_i (mrad A * (Mabs Rm))%Ri i j)
 -> singular A.
Proof.
move=> [Rm [j Hle]].
have Hin: Mabs (mmid A *m (\col_i Rm i j)) <=m: 
  mrad A *m Mabs (\col_i Rm i j).
have He1: ( Mabs (mmid A *m \col_i Rm i j)) =
   (\col_i (Mabs ((mmid A) *  Rm) i j)).
by apply/matrixP => i k;  rewrite !mxE /=; f_equal; 
  apply eq_bigr=> r_; rewrite mxE.
rewrite He1.
have He2: \col_i Mabs (mmid A * Rm) i j =
  \col_i Mabs (1 - (1 - (mmid A) * Rm)) i j.
by apply/matrixP=> i k; rewrite !mxE oppr_sub addrCA addrN addr0.
rewrite He2.
move=> i j0; apply Rle_trans with ( (1 + (Mabs (1 - mmid A * Rm))) i j)%Ri.
rewrite !mxE.  
apply Rle_trans with (1:= Rabs_triang _ _).
rewrite Rabs_right.
apply Rplus_le_compat_l.
by rewrite Rabs_Ropp; right.
case: (i==j);  rewrite /GRing.natmul /=; toR; fourier.
have Hlei := (Hle i ord0).
rewrite mxE  in Hlei.
apply Rle_trans with (1:=Hlei).
rewrite !mxE; right.
by apply eq_bigr => k _; rewrite !mxE .
have Hsa: {vA| inSetm A vA /\ \det vA == 0}.
apply oettli_prager_sing.
exists (\col_i (Rm i j)).
split; [exact Hin|].
apply/matrixP => Hct.
have Hlej:= Hle j ord0.
move: Hlej; rewrite !mxE.
rewrite (eq_refl j).
apply Rlt_not_le.
apply Rle_lt_trans with 0.
right;  apply: big1 => i _ ; rewrite !mxE .
have Hcti := Hct i ord0.
by rewrite !mxE in Hcti ; rewrite Hcti Rabs_R0 mulr0.
apply Rlt_le_trans with 1; first by (toR; fourier).
rewrite -{1}(addr0 1).
apply Rplus_le_compat_l; apply Rabs_pos.
elim Hsa => [va [H1 H2]].
by exists va; split; [|apply/eqP].
Qed.


Corollary cor34:
  \det (mmid A) != 0 -> 1 <=
  \big[Rmax/0]_(j< n) (mrad A * Mabs (mmid A)^-1) j j  
  -> singular A.
Proof.
move=> Hn0 Hle.
apply  thm33.
exists ((mmid A)^-1).
elim: (@max_exists n (Hp n) (fun j =>(mrad A * Mabs (mmid A)^-1) j j)).
move=> j Heq.
exists j.
move=> i k.
rewrite  mulrV.
rewrite addrN.
case C: (i==j).
have C1: i=j.
by move/eqP: C.
rewrite C1 (ord1_0 k) !mxE .
rewrite !mxE in Heq.
rewrite  -Heq.
apply Rle_trans with (2:=Hle).
by right; rewrite Rabs_R0 addr0 eq_refl.
rewrite !mxE C Rabs_R0 addr0.
apply big_rel.
by right.
 move=> _ x _ y Hx Hy.
by rewrite /GRing.natmul -(addr0 0) /=; apply Rplus_le_compat.
move=> k0 _; apply Rmult_pos; rewrite !mxE; [apply rad_pos|apply Rabs_pos].
apply Hn0.
move=> i; rewrite !mxE.
apply big_rel.
by right.
by move=> _ x _ y Hx Hy; apply Rplus_le_le_0_compat.
move=> k _; apply Rmult_pos; rewrite !mxE; [apply rad_pos|apply Rabs_pos].
Qed.


End midInv_RegCrit.


Section eigen_RegCrit.

Variable n': nat.
Notation n:= n'.+1.

Variable A: matrix closedInt n n.

Variable eigens: matrix R n n -> seq R.
Hypothesis HeigAtA: 
forall A lam, lam \in (eigenv A) <-> (lam \in (eigens A)).


Theorem thm41_lam_minmax:
 lam_max eigens ((mrad A) ^T *m mrad A) < 
  lam_min eigens ((mmid A)^T * mmid A) ->
regular A.
Proof.
move=> Hlam.
apply nsing_reg; rewrite sing_ineq_mid_rad => Hsing.
elim Hsing => x [Hle Hdif0].
move/eqP: Hdif0 => Hdif0.
elim: (Rmat_neq Hdif0) => i0 [j0 H].
rewrite mxE in H.
elim (Rlt_not_le _ _ Hlam).

have Hxtx: 0 <=m: (x^T *m x)^-1.
rewrite (mx11_scalar (x^T *m x)) det_scalar_invmx.
move=> i j; rewrite  !mxE (ord1_0 i) (ord1_0 j) /=.
rewrite {1}/GRing.natmul /= /GRing.inv /= /Rinvx.
case C: ((\sum_(k < n) x^T ord0 k * x k ord0) *+ 1 != 0).
by rewrite /GRing.natmul /=; left; apply Rinv_0_lt_compat; apply: sum_sqr_gt.
by rewrite /GRing.natmul /=; left; apply: sum_sqr_gt.
by rewrite -(mx11_scalar); apply vec_prod_unit.

have Hm: @scalar_mx _ 1 (lam_min eigens ((mmid A)^T * mmid A)) <=m: 
    (lam_max eigens ((mrad A)^T *m mrad A))%:M.
have HR := Rayleigh_min Hdif0 (HeigAtA ((mmid A)^T* (mmid A))).
apply Mle_trans with (1:= Rayleigh_min Hdif0 (HeigAtA((mmid A)^T* (mmid A)))).
rewrite mulmxA -trmx_mul.
apply Mle_trans with 
  (Mabs ((mmid A *m x)^T) *m Mabs (mmid A *m x)  / (x^T *m x) ).
apply: Mmul_mle_compat_r; last by done.
apply: Mle_trans (Mabs_mul _ _).

rewrite mulmxA; move=> i j; rewrite !mxE;  apply RRle_abs.

apply Mle_trans with 
  ( (mrad A *m  Mabs  x)^T *m  (mrad A *m  Mabs x) / (x^T *m x) ).
apply: Mmul_mle_compat_r; last by done.
apply: Mmul_mle_compat.
by rewrite /trmx=> i j; move: (Hle j i); rewrite !mxE.
apply Hle.
apply Mabs_pos.
apply Mabs_pos.
rewrite trmx_mul .
apply Mle_trans with 
( ((Mabs x)^T *m ((mrad A)^T *m mrad A) *m Mabs x) / (x^T *m x) ).
by move=> i j; rewrite !mulmxA; right.
have Heqx: x^T *m x = (Mabs x) ^T *m Mabs x.
apply/matrixP=> i j; rewrite !vec_prod_sum !mxE (ord1_0 i) (ord1_0 j) /= !GRing.mulr1n.
apply: eq_bigr => k _; rewrite /trmx !mxE; toR.
have Heqx:  (x k ord0 * x k ord0 = Rsqr (x k ord0))%Re by rewrite //.
have Heqxa:  (Rabs (x k ord0) * Rabs (x k ord0) = Rsqr (Rabs (x k ord0)))%Re 
  by rewrite //.
rewrite Heqx Heqxa; apply Rsqr_abs.
rewrite Heqx.
apply: Rayleigh_max.
apply/eqP=> Habs; move/eqP: Hdif0=> Hdif0 ; 
  apply: Hdif0;  apply/matrixP => i j; rewrite mxE.
move: Habs; rewrite -matrixP=> Habs.
have Hij := Habs i j.
rewrite !mxE in Hij.
move: Hij; rewrite /Rabs ; case C: (Rcase_abs (x i j)) => Hij.
have J: (- - x i j)%Re = (- 0)%Re.
by apply  Ropp_eq_compat.
by rewrite  -(Ropp_involutive (x i j)) J Ropp_0.
done.
done.
have Hm1 := mx11_scalar ((lam_min eigens ((mmid A)^T * mmid A))%:M).
rewrite Hm1  in Hm.
have Hm2 := mx11_scalar ((lam_max eigens ((mrad A)^T * mrad A))%:M).
rewrite Hm2 in Hm.
rewrite !mxE in Hm.
have Hm0:= Hm ord0 ord0.
by rewrite !mxE /= !GRing.mulr1n in Hm0.
Qed.


End eigen_RegCrit.


Section posDef_RegCrit.

Variable n': nat.
Notation n:= n'.+1.

Variable A: matrix closedInt n n.

Variable eigens: matrix R n n -> seq R.
Hypothesis HeigAtA: 
forall A lam, lam \in (eigenv A) <-> (lam \in (eigens A)).

Variable mnorm: can_norm.

Theorem thm51_posdef:
pos_def_gen ((mmid A)^T * mmid A - || (mrad A) ^T * mrad A : mnorm| *m: 1 ) ->
  regular A.
Proof.
move=> Hpos.
apply nsing_reg; rewrite sing_ineq_mid_rad_norm => Hsing.
elim Hsing => x [Hle Hnorm1].
have Hdif0 : x != 0.
apply/vec_prod_pos; rewrite Hnorm1 => i j; 
  rewrite !mxE (ord1_0 i) (ord1_0 j) //= mulr1n; toR; apply Rlt_0_1.
elim: (Rmat_neq Hdif0) => i0 [j0 H].
rewrite mxE in H.
rewrite /pos_def_gen in Hpos.
have Hposx := Hpos x Hdif0.
have Hc:  x^T *m ((mmid A)^T * mmid A - || (mrad A)^T * mrad A : mnorm| *m: 1) *m x
 <=m: 0.
rewrite mulmx_subr mulmx_subl.
have H2 : x^T *m ((mmid A)^T * mmid A) *m x <=m:
  x^T *m (|| (mrad A)^T * mrad A : mnorm| *m: 1) *m x.
rewrite mulmxA -trmx_mul.
apply Mle_trans with 
  (Mabs ((mmid A *m x)^T) *m Mabs (mmid A *m x)  ).
apply: Mle_trans (Mabs_mul _ _).
rewrite mulmxA; move=> i j; rewrite !mxE;  apply RRle_abs.
apply Mle_trans with 
  ( (mrad A *m  Mabs  x)^T *m  (mrad A *m  Mabs x) ).
apply: Mmul_mle_compat.
by rewrite /trmx=> i j; move: (Hle j i); rewrite !mxE.
apply Hle.
apply Mabs_pos.
apply Mabs_pos.
rewrite trmx_mul .
apply Mle_trans with 
( ((Mabs x)^T *m ((mrad A)^T *m mrad A) *m Mabs x) ).
by move=> i j; rewrite !mulmxA; right.
apply Mle_trans with 
  ((Mabs x)^T *m ((mrad A)^T *m mrad A) *m Mabs x / ((Mabs x)^T *m Mabs x)).
by rewrite -vec_prod_abs Hnorm1 invr1 mulr1 => i j; right.
apply Mle_trans with (lam_max eigens (((mrad A)^T *m mrad A)))%:M.
apply: Rayleigh_max.
apply/eqP=> Habs; move/eqP: Hdif0=> Hdif0 ; 
  apply: Hdif0;  apply/matrixP => i j; rewrite mxE.
move: Habs; rewrite -matrixP=> Habs.
have Hij := Habs i j.
rewrite !mxE in Hij.
move: Hij; rewrite /Rabs ; case C: (Rcase_abs (x i j)) => Hij.
have J: (- - x i j)%Re = (- 0)%Re.
by apply  Ropp_eq_compat.
by rewrite  -(Ropp_involutive (x i j)) J Ropp_0.
done.
apply HeigAtA.
apply Mle_trans with (specr eigens ((mrad A)^T *m mrad A))%:M.
move=> i j; rewrite !mxE (ord1_0 i) (ord1_0 j) //= !mulr1n /lam_max /specr.
elim ( (eigens ((mrad A)^T *m mrad A))).
apply: big_rel.
by rewrite  /last Rabs_R0; right.
by apply rmax4.
by move=> *; apply: RRle_abs.
move=> t s IH; rewrite !big_cons.
case Hs: (s != [::]).
apply rmax4.
apply RRle_abs.
by rewrite -(last_cons t Hs).
move/eqP: Hs => ->.
rewrite !big_nil; apply rmax4; apply RRle_abs.

apply Mle_trans with (||(mrad A)^T * mrad A : mnorm|)%:M.
move=> i j; rewrite !mxE  (ord1_0 i) (ord1_0 j) //= !mulr1n; apply specr_le_norm.
by done.
by rewrite -scalemxAr mulmx1 -scalemxAl Hnorm1 scalemx1 => i j; right.
move=> i j; have h3 := H2 i j.
rewrite addrC mxE.
apply Rplus_le_reg_l with (( x^T *m (||(mrad A)^T * mrad A : mnorm| *m: 1) *m x) i j) .
rewrite -Rplus_assoc !mxE; toR; ring_simplify.
move: h3; rewrite !mxE; toR.
by elim (Rlt_not_le _ _ (Hposx ord0 ord0)).
Qed.




End posDef_RegCrit.