with Globals, Strings;
use Globals, Strings;

package body Linear.Contr is

  use SM,SV,SV_Ops,SSV;

  procedure Contr_Matrix(M: in out SMatrix) is
     --- M --> -I-1/(M-I)
  begin
     AddConst(Scal(NegOne),M);
     if Verbosity>1 then Show0("inverting matrix"); end if;
     Invert(M);               -- 1/(DG-I)
     if Verbosity>1 then Show0("done"); end if;
     AddConst(Scal(One),M);   -- 1/(DG-I)+I
     Neg(M);                  -- -I-1/(DG-I)
  end Contr_Matrix;

  procedure Prod(Q: in Modes; M: in SMatrix; F: in Fun; G: in out Fun) is
    --- G := M*F
    U: SVector(Q'Range);
  begin
    Copy(F,G);
    Extract(Q,G,U);
    SetZero(G);
    If not IsZero(U) then
      declare
        V: SVector(Q'Range);
      begin
        Prod(M,U,V);
        AddProd(Q,V,G);
      end;
    end if;
  end Prod;

  procedure ContrTmp0(Q: in Modes; M: in SMatrix; F,Tmp1,Tmp2: in out Fun; Err: out Radius) is
    --- F --> (I+M)Map(F)-MF
    H: Fun renames Tmp1;
    MH: Fun renames Tmp2;
    F0: Fun renames Tmp2;
  begin
    Copy(F,F0);
    Map(F0,F);
    Diff(F,F0,H);       -- F-F0
    Err := MaxNorm(H);
    Prod(Q,M,H,MH);     -- M(F-F0)
    Add(MH,F);          -- F + M(F-F0)
  end ContrTmp0;

  procedure Contr0(Q: in Modes; M: in SMatrix; F: in out Fun; Err: out Radius) is
    --- F --> (I+M)Map(F)-MF
    H,MH: Fun;
    F0: Fun renames MH;
  begin
    Copy(F,F0);
    Map(F0,F);
    Diff(F,F0,H);       -- F-F0
    Err := MaxNorm(H);
    Prod(Q,M,H,MH);     -- M(F-F0)
    Add(MH,F);          -- F + M(F-F0)
  end Contr0;

  procedure Contr(Q: in Modes; M: in SMatrix; F: in Fun; H: in out Fun; Err: out Radius) is
    --- H --> Map(F+(I+M)H)-F-MH
    MH,G: Fun;
  begin
    Prod(Q,M,H,MH);
    Sum(H,MH,G);
    Add(F,G);        -- F+(I+M)H
    Map(G,H);        -- new H
    Sub(H,G);
    Err := MaxNorm(G);
    Sub(F,H);
    Sub(MH,H);       -- H-F-MH
  end Contr;

  procedure DContr(Q: in Modes; M: in SMatrix; H1: in Fun; H2: in out Fun) is
    --- H2 = DMap(I+M)H1-MH1
  begin
    Prod(Q,M,H1,H2);
    if IsZero(H2) then
      DMap(H1,H2);
    else
      declare
        MH,H: Fun;
      begin
        Copy(H2,MH);
        Sum(H1,MH,H);
        DMap(H,H2);
        Sub(MH,H2);
      end;
    end if;
  end DContr;

end Linear.Contr;
