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

use type Ada.Containers.Count_Type;

package Dgemm_Pkg
  with SPARK_Mode => On
is

   type Transpose_Option_Type is (No_Transpose, Transpose, Conjugate_Transpose);

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

   type Matrix_Type is record
      Values                         : Float_Vectors.Vector;
      Leading_Dimension              : Integer;
      Row_Count                      : Integer;
      Column_Count                   : Integer;
   end record;

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

   procedure Double_General_Matrix_Multiply
     (Trans_A : Transpose_Option_Type;
      Trans_B : Transpose_Option_Type;
      M       : Integer;
      N       : Integer;
      K       : Integer;
      Alpha   : Float;
      A       : Matrix_Type;
      B       : Matrix_Type;
      Beta    : Float;
      C       : in out Matrix_Type)
   with
     Pre  =>
       M >= 0
       and then N >= 0
       and then K >= 0
       and then A.Leading_Dimension >=
                  Integer'Max
                    (1,
                     (if Trans_A = No_Transpose then M else K))
       and then B.Leading_Dimension >=
                  Integer'Max
                    (1,
                     (if Trans_B = No_Transpose then K else N))
       and then C.Leading_Dimension >= Integer'Max (1, M)
       and then C.Row_Count >= M
       and then C.Column_Count >= N,
     Post =>
       C.Row_Count = C'Old.Row_Count
       and then C.Column_Count = C'Old.Column_Count
       and then C.Leading_Dimension = C'Old.Leading_Dimension
       and then
         (if M = 0
            or else N = 0
            or else ((Alpha = 0.0 or else K = 0) and then Beta = 1.0)
          then
            C = C'Old)
       and then
         (if M > 0 and then N > 0
            and then Alpha = 0.0
            and then Beta = 0.0
          then
            (for all I in 1 .. M =>
               (for all J in 1 .. N =>
                  Float_Vectors.Element
                    (C.Values,
                     Natural ((I - 1) + (J - 1) * C.Leading_Dimension))
                  = 0.0)))
       and then
         (if M > 0 and then N > 0
            and then Alpha = 0.0
            and then Beta /= 0.0
            and then Beta /= 1.0
          then
            (for all I in 1 .. M =>
               (for all J in 1 .. N =>
                  Float_Vectors.Element
                    (C.Values,
                     Natural ((I - 1) + (J - 1) * C.Leading_Dimension))
                  = Beta
                    * Float_Vectors.Element
                        (C'Old.Values,
                         Natural
                           ((I - 1) + (J - 1) * C'Old.Leading_Dimension)))));

end Dgemm_Pkg;