package body Dgemv_Pkg
  with SPARK_Mode => On
is

   function Logical_X_Index
     (I : Integer; Len : Integer; Inc : Integer) return Natural
   is
   begin
      if Inc > 0 then
         return Natural (I - 1) * Natural (Inc);
      else
         return Natural (Len - I) * Natural (-Inc);
      end if;
   end Logical_X_Index;

   function X_Logical
     (V    : Real_Vector_Type;
      I    : Integer;
      Inc  : Integer;
      Len  : Integer) return Float
   is
      Idx : constant Natural := Logical_X_Index (I, Len, Inc);
   begin
      return Float_Vectors.Element (V.Data, Idx);
   end X_Logical;

   function A_At
     (A   : Real_Matrix_Type;
      I   : Integer;
      J   : Integer;
      Lda : Integer) return Float
   is
      Pos : constant Natural := Natural (I + (J - 1) * Lda - 1);
   begin
      return Float_Vectors.Element (A.Data, Pos);
   end A_At;

   function Partial_Dot
     (Trans : Transpose_Kind_Type;
      A     : Real_Matrix_Type;
      X     : Real_Vector_Type;
      I     : Integer;
      K     : Integer;
      M     : Integer;
      N     : Integer;
      Lda   : Integer;
      Inc_X : Integer) return Float
   is
      pragma Unreferenced (M, N);
      Sum   : Float := 0.0;
      Len_X : constant Integer :=
        (if Trans = No_Transpose then N else M);
   begin
      if K <= 0 then
         return 0.0;
      end if;
      for J in 1 .. K loop
         if Trans = No_Transpose then
            Sum := Sum + A_At (A, I, J, Lda)
                        * X_Logical (X, J, Inc_X, Len_X);
         else
            Sum := Sum + A_At (A, J, I, Lda)
                        * X_Logical (X, J, Inc_X, Len_X);
         end if;
      end loop;
      return Sum;
   end Partial_Dot;

   procedure Dgemv
     (Trans : in     Transpose_Kind_Type;
      M     : in     Integer;
      N     : in     Integer;
      Alpha : in     Float;
      A     : in     Real_Matrix_Type;
      Lda   : in     Integer;
      X     : in     Real_Vector_Type;
      Inc_X : in     Integer;
      Beta  : in     Float;
      Y     : in out Real_Vector_Type;
      Inc_Y : in     Integer)
   is
      Y_Len : constant Integer :=
        (if Trans = No_Transpose then M else N);
      X_Len : constant Integer :=
        (if Trans = No_Transpose then N else M);
   begin
      if M = 0 or else N = 0
        or else (Alpha = 0.0 and then Beta = 1.0)
      then
         return;
      end if;

      for I in 1 .. Y_Len loop
         declare
            Yi : constant Natural :=
              (if Inc_Y > 0 then Natural (I - 1) * Natural (Inc_Y)
               else Natural (Y_Len - I) * Natural (-Inc_Y));
            Y_Old_Val : constant Float :=
              Float_Vectors.Element (Y.Data, Yi);
            Sum     : Float := 0.0;
            New_Val : Float;
         begin
            if Alpha /= 0.0 then
               for J in 1 .. X_Len loop
                  declare
                     Xj : constant Natural :=
                       (if Inc_X > 0 then Natural (J - 1) * Natural (Inc_X)
                        else Natural (X_Len - J) * Natural (-Inc_X));
                     X_Val : constant Float :=
                       Float_Vectors.Element (X.Data, Xj);
                     A_Idx : constant Natural :=
                       (if Trans = No_Transpose then
                           Natural (I + (J - 1) * Lda - 1)
                        else
                           Natural (J + (I - 1) * Lda - 1));
                     A_Val : constant Float :=
                       Float_Vectors.Element (A.Data, A_Idx);
                  begin
                     Sum := Sum + A_Val * X_Val;
                  end;
                  pragma Loop_Invariant (Sum = Sum);
               end loop;
            end if;

            if Beta = 0.0 then
               New_Val := Alpha * Sum;
            else
               New_Val := Beta * Y_Old_Val + Alpha * Sum;
            end if;

            Float_Vectors.Replace_Element (Y.Data, Yi, New_Val);
         end;
      end loop;
   end Dgemv;

end Dgemv_Pkg;