--  Cholesky_Factorization — body for Dpotrf_Pkg.

with Ada.Numerics.Elementary_Functions;

package body Dpotrf_Pkg with SPARK_Mode => On is

   use Ada.Numerics.Elementary_Functions;

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

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

   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)
   is
      Sum_Sq   : Float;
      Sum_Prod : Float;
      Diag_Sq  : Float;
      Diag_Val : Float;
   begin
      --  Argument validation. LAPACK convention: -i for illegal
      --  argument i, in the declared parameter order
      --  (Uplo, N, A, Lda). On illegal argument we must leave A
      --  unchanged (Post: Info.Code < 0 => A = A'Old).

      if Uplo /= UPPER and then Uplo /= LOWER then
         Info.Code := -1;
         return;
      end if;

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

      if Lda < Integer'Max (1, N) then
         Info.Code := -4;
         return;
      end if;

      --  Trivial dimension: nothing to factor; A is left untouched.
      if N = 0 then
         Info.Code := 0;
         return;
      end if;

      --  Column loop — LAPACK Cholesky.

      if Uplo = LOWER then
         for K in 1 .. N loop
            --  Σ_{J=1..K-1} L(K, J)²
            Sum_Sq := 0.0;
            for J in 1 .. K - 1 loop
               declare
                  L_KJ : constant Float := Get_Elem (A, K, J);
               begin
                  Sum_Sq := Sum_Sq + L_KJ * L_KJ;
               end;
               pragma Loop_Invariant (Sum_Sq >= 0.0);
            end loop;

            Diag_Sq := Get_Elem (A, K, K) - Sum_Sq;
            if Diag_Sq <= 0.0 then
               Info.Code := K;
               return;
            end if;

            Diag_Val := Sqrt (Diag_Sq);
            Set_Elem (A, K, K, Diag_Val);

            for I in K + 1 .. N loop
               Sum_Prod := 0.0;
               for J in 1 .. K - 1 loop
                  Sum_Prod := Sum_Prod
                              + Get_Elem (A, I, J) * Get_Elem (A, K, J);
               end loop;
               Set_Elem
                 (A, I, K,
                  (Get_Elem (A, I, K) - Sum_Prod) / Diag_Val);
               pragma Loop_Invariant (Diag_Val > 0.0);
            end loop;
            pragma Loop_Invariant (K <= N);
         end loop;
      else
         for K in 1 .. N loop
            --  Σ_{J=1..K-1} U(J, K)²
            Sum_Sq := 0.0;
            for J in 1 .. K - 1 loop
               declare
                  U_JK : constant Float := Get_Elem (A, J, K);
               begin
                  Sum_Sq := Sum_Sq + U_JK * U_JK;
               end;
               pragma Loop_Invariant (Sum_Sq >= 0.0);
            end loop;

            Diag_Sq := Get_Elem (A, K, K) - Sum_Sq;
            if Diag_Sq <= 0.0 then
               Info.Code := K;
               return;
            end if;

            Diag_Val := Sqrt (Diag_Sq);
            Set_Elem (A, K, K, Diag_Val);

            for J in K + 1 .. N loop
               Sum_Prod := 0.0;
               for I in 1 .. K - 1 loop
                  Sum_Prod := Sum_Prod
                              + Get_Elem (A, I, K) * Get_Elem (A, I, J);
               end loop;
               Set_Elem
                 (A, K, J,
                  (Get_Elem (A, K, J) - Sum_Prod) / Diag_Val);
               pragma Loop_Invariant (Diag_Val > 0.0);
            end loop;
            pragma Loop_Invariant (K <= N);
         end loop;
      end if;

      Info.Code := 0;
   end Compute_Cholesky_Factorization;

end Dpotrf_Pkg;