with Globals, Protected_Counters;
use Globals, Protected_Counters;

pragma Elaborate_All (Globals,Protected_Counters);

package body Vectors.Ops is

  procedure Mult(S: in Scalar; V: in out Vector) is
  begin
    for N in V'First .. EffLast(V) loop
      Mult(S,V(N));
    end loop;
  end Mult;

  procedure Prod(S: in Scalar; V1: in Vector; V2: in out Vector) is
  begin
    if CheckRange and then ((V1'First /= V2'First) or (V1'Last /= V2'Last)) then
      raise Dimension_Error with "Vectors.Prod error";
    end if;
    for N in V2'Range loop
      Prod(S,V1(N),V2(N));
    end loop;
  end Prod;

  function "*"(S: Scalar; V: Vector) return Vector is
    W: Vector(V'Range);
  begin
    for N in V'Range loop
      Prod(S,V(N),W(N));
    end loop;
    return W;
  end "*";

  procedure AddProd(S: in Scalar; V1: in Vector; V2: in out Vector) is
  begin
    if IsZero(S) then return; end if;
    if CheckRange and then ((V1'First /= V2'First) or (V1'Last /= V2'Last)) then
      raise Dimension_Error with "Vectors.AddProd error";
    end if;
    for N in V1'First .. EffLast(V1) loop
      AddProd(S,V1(N),V2(N));
    end loop;
  end AddProd;

  procedure ScalProd(V1,V2: in Vector; S: in out Scalar) is
  begin
    SetZero(S);
    if CheckRange and then ((V1'First /= V2'First) or (V1'Last /= V2'Last)) then
      raise Dimension_Error with "Vectors.ScalProd error";
    end if;
    for N in V2'Range loop
      AddProd(V1(N),V2(N),S);
    end loop;
  end ScalProd;

  function "*"(V1,V2: Vector) return Scalar is
    S: Scalar;
  begin
    ScalProd(V1,V2,S);
    return S;
  end "*";

  procedure Mult(A: in Matrix; V,Tmp: in out Vector) is
  begin
    Copy(V,Tmp);
    Prod(A,Tmp,V);
  end Mult;

  procedure Mult(A: in Matrix; V: in out Vector) is
  begin
    if not IsZero(V) then
      declare
        Tmp: Vector(V'Range);
      begin
        Copy(V,Tmp);
        Prod(A,Tmp,V);
      end;
    end if;
  end Mult;

  procedure Prod(A: in Matrix; V1: in Vector; V2: in out Vector) is
    NT: constant Natural := Reserve_Tasks(V2'Last-V2'First,VProd_Parallel);
    F: constant Integer := V1'First;
    L: constant Integer := V1'Last;
  begin
    if CheckRange and then ((A'First(1) /= V2'First) or (A'Last(1) /= V2'Last)
                              or (A'First(2) /= V1'First) or (A'Last(2) /= V1'Last)) then
      raise Dimension_Error with "Vectors.Prod error";
    end if;
    if NT>0 then
      SetZero(V2);
      AddProd_Parallel(A,V1,V2,NT);
      Free_Tasks(NT);
    else
      for I in V2'Range loop
        declare
          V2I: Component renames V2(I);
        begin
          Prod(A(I,F),V1(F),V2I);
          for J in F+1 .. L loop
            AddProd(A(I,J),V1(J),V2I);
          end loop;
        end;
      end loop;
    end if;
  end Prod;

  procedure TransProd(A: in Matrix; V1: in Vector; V2: in out Vector) is
    F: constant Integer := V1'First;
    L: constant Integer := V1'Last;
  begin
    if CheckRange and then ((A'First(2) /= V2'First) or (A'Last(2) /= V2'Last)
                              or (A'First(1) /= V1'First) or (A'Last(1) /= V1'Last)) then
      raise Dimension_Error with "Vectors.TransProd error";
    end if;
    for I in V2'Range loop
      declare
        V2I: Component renames V2(I);
      begin
        Prod(A(F,I),V1(F),V2I);
        for J in F+1 .. L loop
          AddProd(A(J,I),V1(J),V2I);
        end loop;
      end;
    end loop;
  end TransProd;

  function "*"(A: Matrix; V: Vector) return Vector is
    W: Vector(A'Range(1));
  begin
    Prod(A,V,W);
    return W;
  end "*";

  procedure AddProd(A: in Matrix; V1: in Vector; V2: in out Vector) is
    NT: constant Natural := Reserve_Tasks(V2'Last-V2'First,VProd_Parallel);
  begin
    if CheckRange and then ((A'First(1) /= V2'First) or (A'Last(1) /= V2'Last)
                              or (A'First(2) /= V1'First) or (A'Last(2) /= V1'Last)) then
      raise Dimension_Error with "Vectors.AddProd error";
    end if;
    if NT>0 then
      AddProd_Parallel(A,V1,V2,NT);
      Free_Tasks(NT);
    else
      for I in V2'Range loop
        declare
          V2I: Component renames V2(I);
        begin
          for J in V1'Range loop
            AddProd(A(I,J),V1(J),V2I);
          end loop;
        end;
      end loop;
    end if;
  end AddProd;

  procedure AddProd_Parallel(A: in Matrix; V1: in Vector; V2: in out Vector; NT: in Positive) is
    --- deterministic: tasks working on independent data
    F: constant Integer := V2'First;
    L: constant Integer := V2'Last;
    PC: Protected_Counter(F);

    procedure Prod_Loop is
      I: Integer;
    begin
      Proper_Rounding;
      loop
        PC.Next(I);
        exit when I>L;
        for J in V1'Range loop
          AddProd(A(I,J),V1(J),V2(I));
        end loop;
      end loop;
    exception
      when others => Task_Error := True; raise;
    end Prod_Loop;

    task type Prod_Task_Type is end Prod_Task_Type;
    task body Prod_Task_Type is begin Prod_Loop; end Prod_Task_Type;
  begin
    declare
      Prod_Task: array(1 .. NT) of Prod_Task_Type;
      pragma Warnings (Off,Prod_Task);
    begin
      Prod_Loop;
    end;
    if Task_Error then raise Sorry with "Exception in Vectors.AddProd_Parallel"; end if;
  end AddProd_Parallel;

end Vectors.Ops;
