with Globals, Protected_Counters, Ints;
use Globals, Protected_Counters, Ints;

pragma Elaborate_All (Globals,Protected_Counters);

package body Linear is

   use SM,SV;

   procedure Norm(Q: in Modes; V: in SVector; S: in out Scalar) is
      T,Tmp: Scalar;
   begin
      SetZero(S);
      for I in V'Range loop
         Norm(Q(I),T);
         AddProd(T,V(I),S,Tmp);
      end loop;
   end Norm;

   procedure AddProd(Q: in Modes; V: in SVector; F: in out Fun) is
   begin
      for I in V'Range loop
         AddProd(Q(I),V(I),F,Check_Modes);
      end loop;
   end AddProd;

   procedure Extract(Q: in Modes; F: in out Fun; V: in out SVector) is
   begin
      for I in Q'Range loop
         Extract(Q(I),F,V(I),Check_Modes);
      end loop;
   end Extract;

   procedure Op_Norm(Q: in Modes; S: out Radius) is
      --- deterministic: tasks working on independent data, and PMax is exact
      Small: constant Flt := Flt(7)/Flt(8);
      Task_Error: Boolean := False;
      PC: Protected_Counter(Q'First);

      protected type PMax is procedure Max(T: in Radius); end PMax;
      protected body PMax is procedure Max(T: in Radius) is begin S := RMax(T,S); end Max; end PMax;
      PM: PMax;

      procedure Column_Loop is
         J: Integer;
         S2,T: Radius;
         S1: Scalar;
         F1,F2: Fun;
      begin
         Proper_Rounding;
         loop
            PC.Next(J);
            if (J mod 1000 = 0) and (Verbosity>1) then Show1("j=",J); end if;
            exit when J>Q'Last;
            Assign(Q(J),F1);
            Norm(Q(J),S1);
            LinOp(F1,F2);
            S2 := MaxNorm(F2);
            T := S2/Inf(S1);
            PM.Max(T);
            if Verbosity>2 or else T>Small then
               Show(Q(J));
               Show1("ratio:",T);
            end if;
         end loop;
      exception
         when others => Task_Error := True; raise;
      end Column_Loop;

      task type Column_Task_Type is pragma Storage_Size (Stack_Size); end Column_Task_Type;
      task body Column_Task_Type is begin Column_Loop; end Column_Task_Type;

      NT: constant Natural := Reserve_Tasks(Q'Last-Q'First,Op_Norm_Parallel);
   begin
      S := Zero;
      if NT=0 then
         Column_Loop;
      else
         declare
            Column_Task: array(1 .. NT) of Column_Task_Type;
            pragma Warnings (Off,Column_Task);
         begin
            Column_Loop;
         end;
         Free_Tasks(NT);
         if Task_Error then raise Sorry with "Task_Error"; end if;
      end if;
   end Op_Norm;

   procedure Enclosure(Q: in Modes; A: in out SMatrix) is
      --- deterministic: tasks working on independent data
      Task_Error: Boolean := False;
      PC: Protected_Counter(Q'First);

      procedure Column_Loop is
         J: Integer;
         V: SVector(Q'Range);
         F1,F2: Fun;
      begin
         Proper_Rounding;
         loop
            PC.Next(J);
            exit when J>Q'Last;
            Assign(Q(J),F1);
            LinOp(F1,F2);
            Extract(Q,F2,V);
            for I in Q'Range loop Copy(V(I),A(I,J)); end loop;
            if Verbosity>2 then Show(Q(J)); end if;
         end loop;
      exception
         when others => Task_Error := True; raise;
      end Column_Loop;

      task type Column_Task_Type is pragma Storage_Size (Stack_Size); end Column_Task_Type;
      task body Column_Task_Type is begin Column_Loop; end Column_Task_Type;

      NT: constant Natural := Reserve_Tasks(Q'Last-Q'First,Enclosure_Parallel);
   begin
      if NT=0 then
         Column_Loop;
      else
         declare
            Column_Task: array(1 .. NT) of Column_Task_Type;
            pragma Warnings (Off,Column_Task);
         begin
            Column_Loop;
         end;
         Free_Tasks(NT);
         if Task_Error then raise Sorry with "Task_Error"; end if;
      end if;
   end Enclosure;

   procedure Norm(A: in SMatrix; S: in out Scalar) is
      F1: constant Integer := A'First(1);
      L1: constant Integer := A'Last(1);
      SJ,T: Scalar;
   begin
      SetZero(S);
      for J in reverse A'Range(2) loop
         Norm(A(L1,J),SJ);
         for I in reverse F1 .. L1-1 loop
            Norm(A(I,J),T);
            Add(T,SJ);
         end loop;
         Max(SJ,S);
      end loop;
   end Norm;

   function Norm_Bound(A: SMatrix; B: Scalar) return Boolean is
      --- checks if Norm(A)<=B
      F1: constant Integer := A'First(1);
      L1: constant Integer := A'Last(1);
      SJ,T: Scalar;
   begin
      for J in reverse A'Range(2) loop
         Norm(A(L1,J),SJ);
         for I in reverse F1 .. L1-1 loop
            Norm(A(I,J),T);
            Add(T,SJ);
         end loop;
         if Compare(SJ,B)>0 then return False; end if;
      end loop;
      return True;
   end Norm_Bound;

   procedure OpNorm(Q: in Modes; A: in SMatrix; S: in out Scalar) is
      T,SJ: Scalar;
      W: SVector(Q'Range);
   begin
      SetZero(S);
      for I in Q'Range loop Norm(Q(I),W(I)); end loop;
      for J in reverse A'Range(2) loop
         SetZero(SJ);
         for I in reverse A'Range(1) loop
            Norm(A(I,J),T);
            AddProd(W(I),T,SJ);
         end loop;
         Quot(SJ,W(J),T);
         Max(T,S);
      end loop;
   end OpNorm;

   procedure SpecRadius(A,Tmp: in out SMatrix; S: in out Scalar; Pow2: in Positive := 5) is
      --- A and Tmp get modified
      Show_Norm: constant Boolean := (Verbosity>1);
      T: Scalar;
   begin
      for K in 1 .. Pow2 loop
         Prod(A,A,Tmp);
         Swap(Tmp,A);
         if Show_Norm then
            Norm(A,S);
            Show1("matrix power norm ",S);
         end if;
      end loop;
      if not Show_Norm then Norm(A,S); end if;
      for K in 1 .. Pow2 loop
         Sqrt(S,T);
         Swap(T,S);
      end loop;
   end SpecRadius;

   procedure SpecRadius(A: in SMatrix; S: in out Scalar; Pow2: in Positive := 5) is
      B,Tmp: SMatrix(A'Range(1),A'Range(2));
   begin
      Copy(A,B);
      SpecRadius(B,Tmp,S,Pow2);
   end SpecRadius;

   function SpecRadius_Bound(A: SMatrix; R: Scalar; Pow2: Positive := 5) return Boolean is
      --- check if SpecRadius(A)<=R
      S,T: Scalar;
      B: SMatrix(A'Range(1),A'Range(2));
   begin
      Prod(A,A,B);
      Prod(R,R,S);
      if Norm_Bound(B,S) then return True; end if;
      declare
         Tmp: SMatrix(A'Range(1),A'Range(2));
      begin
         for K in 2 .. Pow2 loop
            Swap(B,Tmp);
            Swap(S,T);
            Prod(Tmp,Tmp,B);
            Prod(T,T,S);
            if Norm_Bound(B,S) then return True; end if;
         end loop;
      end;
      if Verbosity>2 then
         Norm(B,S);
         for K in 1 .. Pow2 loop
            Sqrt(S,T);
            Swap(T,S);
         end loop;
         Show1("radius bound",S);
      end if;
      return False;
   end SpecRadius_Bound;

   procedure Spec_Radius(Q: in Modes; S: in out Scalar; Pow2: in Positive := 5) is
      procedure Enclose is new Enclosure (LinOp => LinOp);
      A,Tmp: SMatrix(Q'Range,Q'Range);
   begin
      Enclose(Q,A);
      SpecRadius(A,Tmp,S,Pow2);
   end Spec_Radius;

   function Spec_Radius_Bound(Q: Modes; R: Scalar; Pow2: Positive := 5) return Boolean is
      --- check if SpecRadius(LinOp)<=R
      procedure Enclose is new Enclosure (LinOp => LinOp);
      A: SMatrix(Q'Range,Q'Range);
   begin
      Enclose(Q,A);
      return SpecRadius_Bound(A,R,Pow2);
   end Spec_Radius_Bound;

end Linear;
