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

use type Ada.Containers.Count_Type;

package Dgetrf_Pkg
  with SPARK_Mode => On
is

   type Info_Code_Type is record
      Value                          : Integer;
   end record;

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

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

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

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

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

   type Pivot_Vector_Type is record
      Item_Count                     : Integer;
      Indices                        : Integer_Vectors.Vector;
   end record;

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

   function Matrix_Element
     (A : Matrix_Type; I, J : Integer) return Float
   is (Float_Vectors.Element
        (A.Elements, Natural ((I - 1) + (J - 1) * A.Leading_Dimension)))
   with Ghost,
        Pre => I in 1 .. A.Row_Count
               and then J in 1 .. A.Column_Count
               and then A.Leading_Dimension >= A.Row_Count
               and then Float_Vectors.Length (A.Elements) >=
                          Ada.Containers.Count_Type
                            (A.Leading_Dimension * A.Column_Count);

   function Pivot_Element
     (IPIV : Pivot_Vector_Type; K : Integer) return Integer
   is (Integer_Vectors.Element (IPIV.Indices, Natural (K - 1)))
   with Ghost,
        Pre => K in 1 .. IPIV.Item_Count
               and then Integer_Vectors.Length (IPIV.Indices) >=
                          Ada.Containers.Count_Type (IPIV.Item_Count);

   function U_Diagonal_Element
     (A : Matrix_Type; K : Integer) return Float
   is (Matrix_Element (A, K, K))
   with Ghost,
        Pre => K in 1 .. Integer'Min (A.Row_Count, A.Column_Count)
               and then A.Leading_Dimension >= A.Row_Count
               and then Float_Vectors.Length (A.Elements) >=
                          Ada.Containers.Count_Type
                            (A.Leading_Dimension * A.Column_Count);

   procedure Compute_Lu_Factorization
     (M    : in     Integer;
      N    : in     Integer;
      A    : in out Matrix_Type;
      LDA  : in     Integer;
      IPIV : in out Pivot_Vector_Type;
      Info :    out Info_Code_Type)
   with
      Pre  => M >= 0
              and then N >= 0
              and then LDA >= Integer'Max (1, M)
              and then A.Row_Count = M
              and then A.Column_Count = N
              and then A.Leading_Dimension = LDA
              and then Float_Vectors.Length (A.Elements) >=
                         Ada.Containers.Count_Type (LDA * N)
              and then Integer_Vectors.Length (IPIV.Indices) >=
                         Ada.Containers.Count_Type (Integer'Min (M, N))
              and then IPIV.Item_Count >= Integer'Min (M, N),
      Post => (if Info.Value < 0 then
                 A = A'Old and then IPIV = IPIV'Old)
              and then
              (if Integer'Min (M, N) = 0 and then Info.Value = 0 then
                 A = A'Old)
              and then
              (if Info.Value >= 0 and then Integer'Min (M, N) > 0 then
                 A.Row_Count = M
                 and then A.Column_Count = N
                 and then A.Leading_Dimension = LDA
                 and then IPIV.Item_Count = Integer'Min (M, N)
                 and then
                   (for all K in 1 .. Integer'Min (M, N) =>
                      Pivot_Element (IPIV, K) >= K
                      and then Pivot_Element (IPIV, K) <= M))
              and then
              (if Info.Value = 0 and then Integer'Min (M, N) > 0 then
                 (for all K in 1 .. Integer'Min (M, N) =>
                    U_Diagonal_Element (A, K) /= 0.0));

end Dgetrf_Pkg;