with SPARK.Containers.Formal.Unbounded_Vectors;
with Ada;
with Ada.Containers;

use type Ada.Containers.Count_Type;

package Dgemv_Pkg
  with SPARK_Mode => On
is

   type Transpose_Kind_Type is (No_Transpose, Transpose, Conjugate_Transpose);

   package Float_Vectors is new SPARK.Containers.Formal.Unbounded_Vectors
      (Index_Type   => Natural,
       Element_Type => Float);

   type Real_Vector_Type is record
      Data                           : Float_Vectors.Vector;
   end record;

   package Real_Vector_Type_Vectors is new SPARK.Containers.Formal.Unbounded_Vectors
      (Index_Type   => Natural,
       Element_Type => Real_Vector_Type);

   type Real_Matrix_Type is record
      Data                           : Float_Vectors.Vector;
      Lda                            : Integer;
   end record;

   package Real_Matrix_Type_Vectors is new SPARK.Containers.Formal.Unbounded_Vectors
      (Index_Type   => Natural,
       Element_Type => Real_Matrix_Type);

   function Logical_X_Index
     (I : Integer; Len : Integer; Inc : Integer) return Natural
   with Ghost,
        Pre  => Inc /= 0 and then Len >= 1 and then I in 1 .. Len,
        Post => Logical_X_Index'Result =
                (if Inc > 0 then Natural (I - 1) * Natural (Inc)
                 else Natural (Len - I) * Natural (-Inc));

   function X_Logical
     (V    : Real_Vector_Type;
      I    : Integer;
      Inc  : Integer;
      Len  : Integer) return Float
   with Ghost,
        Pre => Inc /= 0
               and then Len >= 1
               and then I in 1 .. Len
               and then Float_Vectors.Length (V.Data) >=
                        Ada.Containers.Count_Type
                          (1 + (Len - 1) * (abs Inc));

   function A_At
     (A   : Real_Matrix_Type;
      I   : Integer;
      J   : Integer;
      Lda : Integer) return Float
   with Ghost,
        Pre => I >= 1
               and then J >= 1
               and then Lda >= I
               and then Float_Vectors.Length (A.Data) >=
                        Ada.Containers.Count_Type (I + (J - 1) * Lda);

   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
   with Ghost,
        Pre  => M >= 0
                and then N >= 0
                and then Lda >= Integer'Max (1, M)
                and then Inc_X /= 0
                and then K >= 0
                and then ((Trans = No_Transpose
                           and then I in 1 .. M
                           and then K <= N)
                          or else (Trans /= No_Transpose
                                   and then I in 1 .. N
                                   and then K <= M));

   function Matvec_Element
     (Trans : Transpose_Kind_Type;
      A     : Real_Matrix_Type;
      X     : Real_Vector_Type;
      I     : Integer;
      M     : Integer;
      N     : Integer;
      Lda   : Integer;
      Inc_X : Integer) return Float
   is (Partial_Dot
         (Trans, A, X, I,
          (if Trans = No_Transpose then N else M),
          M, N, Lda, Inc_X))
   with Ghost,
        Pre => M >= 0
               and then N >= 0
               and then Lda >= Integer'Max (1, M)
               and then Inc_X /= 0
               and then ((Trans = No_Transpose and then I in 1 .. M)
                         or else (Trans /= No_Transpose
                                  and then I in 1 .. N));

   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)
   with
      Pre  => M >= 0
              and then N >= 0
              and then Lda >= Integer'Max (1, M)
              and then Inc_X /= 0
              and then Inc_Y /= 0,
      Post =>
        (if M = 0 or else N = 0 or else (Alpha = 0.0 and then Beta = 1.0)
         then Y = Y'Old
         else
           (for all I in 1 ..
              (if Trans = No_Transpose then M else N) =>
                X_Logical
                  (Y, I, Inc_Y,
                   (if Trans = No_Transpose then M else N)) =
                (if Beta = 0.0 then 0.0
                 else Beta * X_Logical
                               (Y'Old, I, Inc_Y,
                                (if Trans = No_Transpose then M else N)))
                + Alpha * Matvec_Element
                            (Trans, A, X, I, M, N, Lda, Inc_X)));

end Dgemv_Pkg;