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

use type Ada.Containers.Count_Type;

package Dgeqrf_Pkg
  with SPARK_Mode => On
is

   type QR_Status_Type is (SUCCESS, ILLEGAL_ARGUMENT_M, ILLEGAL_ARGUMENT_N, ILLEGAL_ARGUMENT_LDA, ILLEGAL_ARGUMENT_LWORK);

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

   type Matrix_Real_Type is record
      Rows                           : Integer;
      Columns                        : Integer;
      Leading_Dimension              : Integer;
      Data                           : Float_Vectors.Vector;
   end record;

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

   type Vector_Real_Type is record
      Item_Count                     : Integer;
      Data                           : Float_Vectors.Vector;
   end record;

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

   function Matrix_Element
     (A : Matrix_Real_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.Rows
            and then J in 1 .. A.Columns
            and then A.Leading_Dimension >= A.Rows
            and then Float_Vectors.Length (A.Data) >=
                       Ada.Containers.Count_Type
                         (A.Leading_Dimension * A.Columns);

   function Tau_Element
     (Tau : Vector_Real_Type; I : Integer) return Float
   is (Float_Vectors.Element (Tau.Data, Natural (I - 1)))
   with Ghost,
     Pre => I in 1 .. Tau.Item_Count
            and then Float_Vectors.Length (Tau.Data) >=
                       Ada.Containers.Count_Type (Tau.Item_Count);

   function Reflector_V_Component
     (A : Matrix_Real_Type; I, K, M_Rows : Integer) return Float
   is (if K < I then 0.0
    elsif K = I then 1.0
    elsif K <= M_Rows then Matrix_Element (A, K, I)
    else 0.0)
   with Ghost,
     Pre => I in 1 .. Integer'Min (A.Rows, A.Columns)
            and then K in 1 .. M_Rows
            and then M_Rows = A.Rows
            and then A.Leading_Dimension >= A.Rows
            and then Float_Vectors.Length (A.Data) >=
                       Ada.Containers.Count_Type
                         (A.Leading_Dimension * A.Columns);

   function Sum_Of_Reflector_Squares
     (A : Matrix_Real_Type; I, K, M_Rows : Integer) return Float
   is (if K > M_Rows then 0.0
    else Reflector_V_Component (A, I, K, M_Rows)
         * Reflector_V_Component (A, I, K, M_Rows)
         + Sum_Of_Reflector_Squares (A, I, K + 1, M_Rows))
   with Ghost,
     Subprogram_Variant => (Decreases => M_Rows - K + 1),
     Pre => I in 1 .. Integer'Min (A.Rows, A.Columns)
            and then K in 1 .. M_Rows + 1
            and then M_Rows = A.Rows
            and then A.Leading_Dimension >= A.Rows
            and then Float_Vectors.Length (A.Data) >=
                       Ada.Containers.Count_Type
                         (A.Leading_Dimension * A.Columns);

   function Reflector_Norm_Squared
     (A : Matrix_Real_Type; I, M_Rows : Integer) return Float
   is (Sum_Of_Reflector_Squares (A, I, I, M_Rows))
   with Ghost;

   function Encodes_Orthogonal_Reflector
     (A : Matrix_Real_Type; Tau : Vector_Real_Type;
   I, M_Rows : Integer) return Boolean
   is (Tau_Element (Tau, I) = 0.0
    or else abs (Tau_Element (Tau, I)
                 * Reflector_Norm_Squared (A, I, M_Rows)
                 - 2.0) < 1.0e-4)
   with Ghost;

   function R_Factor_Element
     (A     : Matrix_Real_Type;
   I, J  : Integer;
   K_Min : Integer) return Float
   is (if I > K_Min then 0.0
    elsif I > J then 0.0
    else Matrix_Element (A, I, J))
   with Ghost;

   procedure Compute_Qr_Factorization
     (A      : in out Matrix_Real_Type;
   Tau    : in out Vector_Real_Type;
   Work   : in out Vector_Real_Type;
   M_Rows : in     Integer;
   N_Cols : in     Integer;
   LDA    : in     Integer;
   LWork  : in     Integer;
   Status :    out QR_Status_Type)
   with
     Pre  => M_Rows >= 0
          and then N_Cols >= 0
          and then LDA >= Integer'Max (1, M_Rows)
          and then (LWork = -1
                    or else (Integer'Min (M_Rows, N_Cols) = 0
                             and then LWork >= 1)
                    or else (M_Rows >= 1 and then LWork >= N_Cols)),
     Post => (if LWork = -1 then
             Status = SUCCESS
             and then A = A'Old
             and then Tau = Tau'Old)
          and then
          (if Status /= SUCCESS then
             A = A'Old and then Tau = Tau'Old)
          and then
          (if Status = SUCCESS
              and then LWork /= -1
              and then Integer'Min (M_Rows, N_Cols) = 0
           then
             A = A'Old and then Tau = Tau'Old)
          and then
          (if Status = SUCCESS
              and then LWork /= -1
              and then Integer'Min (M_Rows, N_Cols) > 0
           then
             A.Rows = M_Rows
             and then A.Columns = N_Cols
             and then A.Leading_Dimension = LDA
             and then Tau.Item_Count = Integer'Min (M_Rows, N_Cols)
             and then
               (for all I in 1 .. Integer'Min (M_Rows, N_Cols) =>
                  Reflector_V_Component (A, I, I, M_Rows) = 1.0)
             and then
               (for all I in 1 .. Integer'Min (M_Rows, N_Cols) =>
                  (for all K in 1 .. I - 1 =>
                     Reflector_V_Component (A, I, K, M_Rows) = 0.0))
             and then
               (for all I in 2 .. Integer'Min (M_Rows, N_Cols) =>
                  (for all J in 1 .. I - 1 =>
                     R_Factor_Element
                       (A, I, J, Integer'Min (M_Rows, N_Cols)) = 0.0))
             and then
               (for all I in 1 .. Integer'Min (M_Rows, N_Cols) =>
                  Encodes_Orthogonal_Reflector (A, Tau, I, M_Rows)));

end Dgeqrf_Pkg;