package body Dgemm_Pkg
  with SPARK_Mode => On
is

   procedure Double_General_Matrix_Multiply
     (Trans_A : Transpose_Option_Type;
      Trans_B : Transpose_Option_Type;
      M       : Integer;
      N       : Integer;
      K       : Integer;
      Alpha   : Float;
      A       : Matrix_Type;
      B       : Matrix_Type;
      Beta    : Float;
      C       : in out Matrix_Type)
   is
      A_Ld : constant Integer := A.Leading_Dimension;
      B_Ld : constant Integer := B.Leading_Dimension;
      C_Ld : constant Integer := C.Leading_Dimension;
   begin
      --  Quick return when either output dimension is empty.
      if M = 0 or else N = 0 then
         return;
      end if;

      --  Quick return when the result equals C unchanged.
      if (Alpha = 0.0 or else K = 0) and then Beta = 1.0 then
         return;
      end if;

      --  Alpha = 0 case: C := Beta * C (or zero).  This branch
      --  discharges the two quantified Post clauses for Alpha = 0.
      if Alpha = 0.0 then
         if Beta = 0.0 then
            for J in 1 .. N loop
               pragma Loop_Invariant
                 (C.Row_Count = C'Loop_Entry.Row_Count);
               pragma Loop_Invariant
                 (C.Column_Count = C'Loop_Entry.Column_Count);
               pragma Loop_Invariant
                 (C.Leading_Dimension = C'Loop_Entry.Leading_Dimension);
               for I in 1 .. M loop
                  pragma Loop_Invariant
                    (C.Row_Count = C'Loop_Entry.Row_Count);
                  pragma Loop_Invariant
                    (C.Column_Count = C'Loop_Entry.Column_Count);
                  pragma Loop_Invariant
                    (C.Leading_Dimension = C'Loop_Entry.Leading_Dimension);
                  Float_Vectors.Replace_Element
                    (C.Values,
                     Natural ((I - 1) + (J - 1) * C_Ld),
                     0.0);
               end loop;
            end loop;
         else
            for J in 1 .. N loop
               pragma Loop_Invariant
                 (C.Row_Count = C'Loop_Entry.Row_Count);
               pragma Loop_Invariant
                 (C.Column_Count = C'Loop_Entry.Column_Count);
               pragma Loop_Invariant
                 (C.Leading_Dimension = C'Loop_Entry.Leading_Dimension);
               for I in 1 .. M loop
                  pragma Loop_Invariant
                    (C.Row_Count = C'Loop_Entry.Row_Count);
                  pragma Loop_Invariant
                    (C.Column_Count = C'Loop_Entry.Column_Count);
                  pragma Loop_Invariant
                    (C.Leading_Dimension = C'Loop_Entry.Leading_Dimension);
                  declare
                     Idx     : constant Natural :=
                       Natural ((I - 1) + (J - 1) * C_Ld);
                     Old_Val : constant Float :=
                       Float_Vectors.Element (C.Values, Idx);
                  begin
                     Float_Vectors.Replace_Element
                       (C.Values, Idx, Beta * Old_Val);
                  end;
               end loop;
            end loop;
         end if;
         return;
      end if;

      --  General case (Alpha /= 0): first scale C by Beta.
      if Beta = 0.0 then
         for J in 1 .. N loop
            pragma Loop_Invariant
              (C.Row_Count = C'Loop_Entry.Row_Count);
            pragma Loop_Invariant
              (C.Column_Count = C'Loop_Entry.Column_Count);
            pragma Loop_Invariant
              (C.Leading_Dimension = C'Loop_Entry.Leading_Dimension);
            for I in 1 .. M loop
               pragma Loop_Invariant
                 (C.Row_Count = C'Loop_Entry.Row_Count);
               pragma Loop_Invariant
                 (C.Column_Count = C'Loop_Entry.Column_Count);
               pragma Loop_Invariant
                 (C.Leading_Dimension = C'Loop_Entry.Leading_Dimension);
               Float_Vectors.Replace_Element
                 (C.Values,
                  Natural ((I - 1) + (J - 1) * C_Ld),
                  0.0);
            end loop;
         end loop;
      elsif Beta /= 1.0 then
         for J in 1 .. N loop
            pragma Loop_Invariant
              (C.Row_Count = C'Loop_Entry.Row_Count);
            pragma Loop_Invariant
              (C.Column_Count = C'Loop_Entry.Column_Count);
            pragma Loop_Invariant
              (C.Leading_Dimension = C'Loop_Entry.Leading_Dimension);
            for I in 1 .. M loop
               pragma Loop_Invariant
                 (C.Row_Count = C'Loop_Entry.Row_Count);
               pragma Loop_Invariant
                 (C.Column_Count = C'Loop_Entry.Column_Count);
               pragma Loop_Invariant
                 (C.Leading_Dimension = C'Loop_Entry.Leading_Dimension);
               declare
                  Idx     : constant Natural :=
                    Natural ((I - 1) + (J - 1) * C_Ld);
                  Old_Val : constant Float :=
                    Float_Vectors.Element (C.Values, Idx);
               begin
                  Float_Vectors.Replace_Element
                    (C.Values, Idx, Beta * Old_Val);
               end;
            end loop;
         end loop;
      end if;

      if K = 0 then
         return;
      end if;

      --  Add alpha * op(A) * op(B) to C.  For real Float, the
      --  Conjugate_Transpose case is identical to Transpose.
      case Trans_A is
         when No_Transpose =>
            case Trans_B is
               when No_Transpose =>
                  --  C := C + alpha * A * B
                  for J in 1 .. N loop
                     pragma Loop_Invariant
                       (C.Row_Count = C'Loop_Entry.Row_Count);
                     pragma Loop_Invariant
                       (C.Column_Count = C'Loop_Entry.Column_Count);
                     pragma Loop_Invariant
                       (C.Leading_Dimension =
                          C'Loop_Entry.Leading_Dimension);
                     for L in 1 .. K loop
                        pragma Loop_Invariant
                          (C.Row_Count = C'Loop_Entry.Row_Count);
                        pragma Loop_Invariant
                          (C.Column_Count = C'Loop_Entry.Column_Count);
                        pragma Loop_Invariant
                          (C.Leading_Dimension =
                             C'Loop_Entry.Leading_Dimension);
                        declare
                           B_Val : constant Float :=
                             Float_Vectors.Element
                               (B.Values,
                                Natural ((L - 1) + (J - 1) * B_Ld));
                           Temp  : constant Float := Alpha * B_Val;
                        begin
                           for I in 1 .. M loop
                              pragma Loop_Invariant
                                (C.Row_Count = C'Loop_Entry.Row_Count);
                              pragma Loop_Invariant
                                (C.Column_Count =
                                   C'Loop_Entry.Column_Count);
                              pragma Loop_Invariant
                                (C.Leading_Dimension =
                                   C'Loop_Entry.Leading_Dimension);
                              declare
                                 C_Idx : constant Natural :=
                                   Natural ((I - 1) + (J - 1) * C_Ld);
                                 A_Val : constant Float :=
                                   Float_Vectors.Element
                                     (A.Values,
                                      Natural
                                        ((I - 1) + (L - 1) * A_Ld));
                                 C_Val : constant Float :=
                                   Float_Vectors.Element
                                     (C.Values, C_Idx);
                              begin
                                 Float_Vectors.Replace_Element
                                   (C.Values,
                                    C_Idx,
                                    C_Val + Temp * A_Val);
                              end;
                           end loop;
                        end;
                     end loop;
                  end loop;

               when others =>
                  --  C := C + alpha * A * op(B), op = transpose
                  for J in 1 .. N loop
                     pragma Loop_Invariant
                       (C.Row_Count = C'Loop_Entry.Row_Count);
                     pragma Loop_Invariant
                       (C.Column_Count = C'Loop_Entry.Column_Count);
                     pragma Loop_Invariant
                       (C.Leading_Dimension =
                          C'Loop_Entry.Leading_Dimension);
                     for L in 1 .. K loop
                        pragma Loop_Invariant
                          (C.Row_Count = C'Loop_Entry.Row_Count);
                        pragma Loop_Invariant
                          (C.Column_Count = C'Loop_Entry.Column_Count);
                        pragma Loop_Invariant
                          (C.Leading_Dimension =
                             C'Loop_Entry.Leading_Dimension);
                        declare
                           B_Val : constant Float :=
                             Float_Vectors.Element
                               (B.Values,
                                Natural ((J - 1) + (L - 1) * B_Ld));
                           Temp  : constant Float := Alpha * B_Val;
                        begin
                           for I in 1 .. M loop
                              pragma Loop_Invariant
                                (C.Row_Count = C'Loop_Entry.Row_Count);
                              pragma Loop_Invariant
                                (C.Column_Count =
                                   C'Loop_Entry.Column_Count);
                              pragma Loop_Invariant
                                (C.Leading_Dimension =
                                   C'Loop_Entry.Leading_Dimension);
                              declare
                                 C_Idx : constant Natural :=
                                   Natural ((I - 1) + (J - 1) * C_Ld);
                                 A_Val : constant Float :=
                                   Float_Vectors.Element
                                     (A.Values,
                                      Natural
                                        ((I - 1) + (L - 1) * A_Ld));
                                 C_Val : constant Float :=
                                   Float_Vectors.Element
                                     (C.Values, C_Idx);
                              begin
                                 Float_Vectors.Replace_Element
                                   (C.Values,
                                    C_Idx,
                                    C_Val + Temp * A_Val);
                              end;
                           end loop;
                        end;
                     end loop;
                  end loop;
            end case;

         when others =>
            --  Trans_A is Transpose or Conjugate_Transpose.
            case Trans_B is
               when No_Transpose =>
                  --  C := C + alpha * A^T * B
                  for J in 1 .. N loop
                     pragma Loop_Invariant
                       (C.Row_Count = C'Loop_Entry.Row_Count);
                     pragma Loop_Invariant
                       (C.Column_Count = C'Loop_Entry.Column_Count);
                     pragma Loop_Invariant
                       (C.Leading_Dimension =
                          C'Loop_Entry.Leading_Dimension);
                     for I in 1 .. M loop
                        pragma Loop_Invariant
                          (C.Row_Count = C'Loop_Entry.Row_Count);
                        pragma Loop_Invariant
                          (C.Column_Count = C'Loop_Entry.Column_Count);
                        pragma Loop_Invariant
                          (C.Leading_Dimension =
                             C'Loop_Entry.Leading_Dimension);
                        declare
                           Temp : Float := 0.0;
                        begin
                           for L in 1 .. K loop
                              declare
                                 A_Val : constant Float :=
                                   Float_Vectors.Element
                                     (A.Values,
                                      Natural
                                        ((L - 1) + (I - 1) * A_Ld));
                                 B_Val : constant Float :=
                                   Float_Vectors.Element
                                     (B.Values,
                                      Natural
                                        ((L - 1) + (J - 1) * B_Ld));
                              begin
                                 Temp := Temp + A_Val * B_Val;
                              end;
                           end loop;
                           declare
                              C_Idx : constant Natural :=
                                Natural ((I - 1) + (J - 1) * C_Ld);
                              C_Val : constant Float :=
                                Float_Vectors.Element (C.Values, C_Idx);
                           begin
                              Float_Vectors.Replace_Element
                                (C.Values,
                                 C_Idx,
                                 C_Val + Alpha * Temp);
                           end;
                        end;
                     end loop;
                  end loop;

               when others =>
                  --  C := C + alpha * A^T * B^T
                  for J in 1 .. N loop
                     pragma Loop_Invariant
                       (C.Row_Count = C'Loop_Entry.Row_Count);
                     pragma Loop_Invariant
                       (C.Column_Count = C'Loop_Entry.Column_Count);
                     pragma Loop_Invariant
                       (C.Leading_Dimension =
                          C'Loop_Entry.Leading_Dimension);
                     for I in 1 .. M loop
                        pragma Loop_Invariant
                          (C.Row_Count = C'Loop_Entry.Row_Count);
                        pragma Loop_Invariant
                          (C.Column_Count = C'Loop_Entry.Column_Count);
                        pragma Loop_Invariant
                          (C.Leading_Dimension =
                             C'Loop_Entry.Leading_Dimension);
                        declare
                           Temp : Float := 0.0;
                        begin
                           for L in 1 .. K loop
                              declare
                                 A_Val : constant Float :=
                                   Float_Vectors.Element
                                     (A.Values,
                                      Natural
                                        ((L - 1) + (I - 1) * A_Ld));
                                 B_Val : constant Float :=
                                   Float_Vectors.Element
                                     (B.Values,
                                      Natural
                                        ((J - 1) + (L - 1) * B_Ld));
                              begin
                                 Temp := Temp + A_Val * B_Val;
                              end;
                           end loop;
                           declare
                              C_Idx : constant Natural :=
                                Natural ((I - 1) + (J - 1) * C_Ld);
                              C_Val : constant Float :=
                                Float_Vectors.Element (C.Values, C_Idx);
                           begin
                              Float_Vectors.Replace_Element
                                (C.Values,
                                 C_Idx,
                                 C_Val + Alpha * Temp);
                           end;
                        end;
                     end loop;
                  end loop;
            end case;
      end case;
   end Double_General_Matrix_Multiply;

end Dgemm_Pkg;