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

use type Ada.Containers.Count_Type;

package Dpotrf_Pkg
  with SPARK_Mode => On
is

   type Triangle_Selector_Type is (UPPER, LOWER);

   type Factorization_Status_Type is record
      Code                           : Integer;
   end record;

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

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

   type Matrix_Storage_Type is record
      Leading_Dimension              : Integer;
      Order                          : Integer;
      Data                           : Float_Vectors.Vector;
   end record;

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

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

   function Lower_Cholesky_Sum_Of_Products
     (A : Matrix_Storage_Type; I, J, K, J_End : Integer) return Float
   is (if K > J_End then 0.0
    else Matrix_Element (A, I, K) * Matrix_Element (A, J, K)
         + Lower_Cholesky_Sum_Of_Products (A, I, J, K + 1, J_End))
   with Ghost,
     Subprogram_Variant => (Decreases => J_End - K + 1),
     Pre => I in 1 .. A.Order
            and then J in 1 .. A.Order
            and then K in 1 .. J_End + 1
            and then J_End in 0 .. Integer'Min (I, J)
            and then A.Leading_Dimension >= A.Order
            and then Float_Vectors.Length (A.Data) >=
                       Ada.Containers.Count_Type
                         (A.Leading_Dimension * A.Order);

   function Upper_Cholesky_Sum_Of_Products
     (A : Matrix_Storage_Type; I, J, K, I_End : Integer) return Float
   is (if K > I_End then 0.0
    else Matrix_Element (A, K, I) * Matrix_Element (A, K, J)
         + Upper_Cholesky_Sum_Of_Products (A, I, J, K + 1, I_End))
   with Ghost,
     Subprogram_Variant => (Decreases => I_End - K + 1),
     Pre => I in 1 .. A.Order
            and then J in 1 .. A.Order
            and then K in 1 .. I_End + 1
            and then I_End in 0 .. Integer'Min (I, J)
            and then A.Leading_Dimension >= A.Order
            and then Float_Vectors.Length (A.Data) >=
                       Ada.Containers.Count_Type
                         (A.Leading_Dimension * A.Order);

   function Lower_Cholesky_Reconstructs
     (A_New : Matrix_Storage_Type;
   A_Old : Matrix_Storage_Type;
   I, J  : Integer)
   return Boolean
   is (abs (Matrix_Element (A_Old, I, J)
         - Lower_Cholesky_Sum_Of_Products (A_New, I, J, 1, J))
    < 1.0e-4)
   with Ghost,
     Pre => I in 1 .. A_New.Order
            and then J in 1 .. I
            and then A_New.Order = A_Old.Order
            and then A_New.Leading_Dimension = A_Old.Leading_Dimension
            and then A_New.Leading_Dimension >= A_New.Order
            and then Float_Vectors.Length (A_New.Data) >=
                       Ada.Containers.Count_Type
                         (A_New.Leading_Dimension * A_New.Order)
            and then Float_Vectors.Length (A_Old.Data) >=
                       Ada.Containers.Count_Type
                         (A_Old.Leading_Dimension * A_Old.Order);

   function Upper_Cholesky_Reconstructs
     (A_New : Matrix_Storage_Type;
   A_Old : Matrix_Storage_Type;
   I, J  : Integer)
   return Boolean
   is (abs (Matrix_Element (A_Old, I, J)
         - Upper_Cholesky_Sum_Of_Products (A_New, I, J, 1, I))
    < 1.0e-4)
   with Ghost,
     Pre => I in 1 .. A_New.Order
            and then J in I .. A_New.Order
            and then A_New.Order = A_Old.Order
            and then A_New.Leading_Dimension = A_Old.Leading_Dimension
            and then A_New.Leading_Dimension >= A_New.Order
            and then Float_Vectors.Length (A_New.Data) >=
                       Ada.Containers.Count_Type
                         (A_New.Leading_Dimension * A_New.Order)
            and then Float_Vectors.Length (A_Old.Data) >=
                       Ada.Containers.Count_Type
                         (A_Old.Leading_Dimension * A_Old.Order);

   procedure Compute_Cholesky_Factorization
     (Uplo : in     Triangle_Selector_Type;
   N    : in     Integer;
   A    : in out Matrix_Storage_Type;
   Lda  : in     Integer;
   Info :    out Factorization_Status_Type)
   with
     Pre  => N >= 0
          and then Lda >= Integer'Max (1, N)
          and then A.Order = N
          and then A.Leading_Dimension = Lda
          and then Float_Vectors.Length (A.Data) >=
                     Ada.Containers.Count_Type
                       (Integer'Max (1, Lda) * Integer'Max (1, N)),
     Post => (if Info.Code < 0 then A = A'Old)
          and then
          (if N = 0 and then Info.Code = 0 then A = A'Old)
          and then
          (if Info.Code = 0 and then N > 0 then
             A.Order = N
             and then A.Leading_Dimension = Lda
             and then
               (if Uplo = LOWER then
                  (for all I in 1 .. N =>
                     (for all J in 1 .. I =>
                        Lower_Cholesky_Reconstructs (A, A'Old, I, J))))
             and then
               (if Uplo = UPPER then
                  (for all I in 1 .. N =>
                     (for all J in I .. N =>
                        Upper_Cholesky_Reconstructs (A, A'Old, I, J)))));

end Dpotrf_Pkg;