package body Dgetrf_Pkg with SPARK_Mode => On is

   function Get_Elem
     (A : Matrix_Type; I, J : Integer) return Float
   is (Float_Vectors.Element
        (A.Elements, Natural ((I - 1) + (J - 1) * A.Leading_Dimension)));

   procedure Set_Elem
     (A : in out Matrix_Type; I, J : Integer; Val : Float) is
   begin
      Float_Vectors.Replace_Element
        (A.Elements, Natural ((I - 1) + (J - 1) * A.Leading_Dimension), Val);
   end Set_Elem;

   procedure Set_Pivot
     (IPIV : in out Pivot_Vector_Type; K, Val : Integer) is
   begin
      Integer_Vectors.Replace_Element
        (IPIV.Indices, Natural (K - 1), Val);
   end Set_Pivot;

   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)
   is
      K_Min   : constant Integer := Integer'Min (M, N);
      P_Idx   : Integer;
      Max_Val : Float;
      Cur     : Float;
      Pivot   : Float;
      Tmp     : Float;
      Recorded_Failure : Boolean := False;
   begin

      if M < 0 then
         Info.Value := -1;
         return;
      end if;

      if N < 0 then
         Info.Value := -2;
         return;
      end if;

      if LDA < Integer'Max (1, M) then
         Info.Value := -4;
         return;
      end if;

      if K_Min = 0 then
         Info.Value := 0;
         return;
      end if;

      A.Row_Count         := M;
      A.Column_Count      := N;
      A.Leading_Dimension := LDA;

      while Natural (Integer_Vectors.Length (IPIV.Indices)) < K_Min loop
         Integer_Vectors.Append (IPIV.Indices, 0);
      end loop;
      IPIV.Item_Count := K_Min;

      Info.Value := 0;

      for K in 1 .. K_Min loop

         P_Idx   := K;
         Max_Val := abs (Get_Elem (A, K, K));
         for I in K + 1 .. M loop
            Cur := abs (Get_Elem (A, I, K));
            if Cur > Max_Val then
               P_Idx   := I;
               Max_Val := Cur;
            end if;
            pragma Loop_Invariant (P_Idx in K .. M);
         end loop;

         Set_Pivot (IPIV, K, P_Idx);

         if P_Idx /= K then
            for J in 1 .. N loop
               Tmp := Get_Elem (A, K, J);
               Set_Elem (A, K, J, Get_Elem (A, P_Idx, J));
               Set_Elem (A, P_Idx, J, Tmp);
            end loop;
         end if;

         Pivot := Get_Elem (A, K, K);
         if Pivot = 0.0 then
            if not Recorded_Failure then
               Info.Value := K;
               Recorded_Failure := True;
            end if;
         else
            for I in K + 1 .. M loop
               Set_Elem (A, I, K, Get_Elem (A, I, K) / Pivot);
            end loop;

            for J in K + 1 .. N loop
               declare
                  U_KJ : constant Float := Get_Elem (A, K, J);
               begin
                  for I in K + 1 .. M loop
                     Set_Elem
                       (A, I, J,
                        Get_Elem (A, I, J)
                        - Get_Elem (A, I, K) * U_KJ);
                  end loop;
               end;
            end loop;
         end if;

         pragma Loop_Invariant
           (IPIV.Item_Count = K_Min
            and then A.Row_Count = M
            and then A.Column_Count = N
            and then A.Leading_Dimension = LDA);
      end loop;

   end Compute_Lu_Factorization;

end Dgetrf_Pkg;