Re: C to Ada conversion

["Dr. Michael Paus" <paus@ib-paus.com>]

tmoran@acm.org wrote:
>>>ok but otherwise it is very dangerous because it does not take into
>>>account that in Ada the index range of arrays does not necessarly
> 
>   On second thought, a more accurate translation from the C, and
> one that would be less general and thus perhaps show better in a
> benchmark, would be to replace
>   type Matrices is array (Integer range <>, Integer range <>) of Float;
> by
>   Rs : constant := 6;
>   type Matrices is array (1 .. Rs, 1 .. Rs) of Float;

Do you think that's fair ?-) The original profile of the C routine was

int MatrixInversion(float* A, int n, float* AInverse)
   // A = input matrix (n x n)
   // n = dimension of A
   // AInverse = inverted matrix (n x n)
   // This function inverts a matrix based on the Gauss Jordan method.
   // The function returns 1 on success, 0 on failure.

so it can in principle be used for any size of array as long as the
dimension n is used consistently for both arrays. There is no implicit
assumption on the size of n (should be > 0 :-) ) and no assumption about
index ranges (where you do not have any choice in C anyway). So, in a
fair comparison the Ada procedure should not make any assumptions which
are more restrictive than the ones in C. Don't you think so too?

Also in your code you write:

   if Det = 0.0 then          -- A pointless test in Ada since an
     raise Constraint_Error;  -- attempted divide by zero will raise
   end if;                    -- Constraint_Error anyway.

This is NOT pointless because a piece of code should also
work correctly if someone turns off the automatic checks in Ada.
The program logic should never rely on these automatic tests.
Your further comment that a test for exactly 0 is not very usefull
is right of course.

When it comes to performance it is important to eliminate as many
constraint checks as possible (without using the brute force method
of just switching them off of course). What do you think about my
quick and dirty proposal? The problems with the index ranges and
the index checks can be solved by introducing an internal method
which maps the index ranges to one well defined range and thus
makes all index checks unneccessary because the compiler (if it
is clever) can statically prove that all loop indices will always
stay in their valid range and there is thus no need for any checks.
In this case you could even switch off all constraint checks
because, if I haven't missed something, all necessary checks are
done explicitly. Also making the internal index range start with
0 instead of 1 might help depending how the compiler computes the
offset to get a certain element in an array.

procedure Inv2 is

    type Matrices is array (Integer range <>, Integer range <>) of Float;

    procedure Matrix_Inversion (
          A       : in out Matrices;
          Inverse :    out Matrices  ) is

       subtype Internal_Index_Type is Integer range 0 .. A'Length(1) - 1;

       subtype Internal_Matrices is Matrices(Internal_Index_Type, Internal_Index_Type);

       procedure Internal_Matrix_Inversion (
             A       : in out Internal_Matrices;
             Inverse :    out Internal_Matrices  ) is
          -- This function inverts a matrix based on the Gauss Jordan method.
          -- It raises Constraint_Error if A or Inverse is not square
          -- or if their index ranges don't match
          -- or if the determinant of A is too small.
          -- Note that input matrix A will be modified.
          Det    : Float := 1.0;
          Temp,
          Pivot,
          Factor : Float;
          -- Ac     : Matrices := A;   -- Note: never used, but in C code
       begin

          Inverse := (others => (others => 0.0));
          for I in Internal_Index_Type loop
             Inverse(I, I) := 1.0;
          end loop;

          for Ipass in Internal_Index_Type loop
             -- For each pass, first find maximum element in the pivot column.
             declare
                Max_Row : Internal_Index_Type := Ipass;
             begin
                for Row in Ipass .. Internal_Index_Type'Last loop
                   if abs(A(Row, Ipass)) > abs(A(Max_Row, Ipass)) then
                      Max_Row := Row;
                   end if;
                end loop;
                -- Interchange the elements of row iPass and Max_Row in both A and Inverse.
                if Max_Row /= Ipass then
                   for Col in Internal_Index_Type loop
                      Temp := Inverse(Ipass, Col);
                      Inverse(Ipass, Col) := Inverse(Max_Row, Col);
                      Inverse(Max_Row, Col) := Temp;
                   end loop;
                end if;
             end;

             -- The current pivot is now A(iPass, iPass).
             -- The determinant is the product of the pivot elements.
             Pivot := A(Ipass, Ipass);
             Det := Det * Pivot;

             if Det = 0.0 then          -- A pointless test in Ada since an
                raise Constraint_Error;  -- attempted divide by zero will raise
             end if;                    -- Constraint_Error anyway.
             -- An erroneous test in C since a very
             -- small, but non-zero, Det will not be
             -- caught, but will cause divide by zero.

             -- Normalize the pivot row by dividing by the pivot element.
             for Col in Internal_Index_Type loop
                Inverse(Ipass, Col) := Inverse(Ipass, Col) / Pivot;
                if Col >= Ipass then
                   A(Ipass, Col) := A(Ipass, Col) / Pivot;
                end if;
             end loop;

             -- Add a multiple of the pivot row to each row.  The multiple factor
             -- is chosen so that the element of A on the pivot column is 0.
             for Row in Internal_Index_Type loop
                -- Note: this loop had two "if (irow != iPass)" statements,
                -- one inside and one outside the column loop.  I've combined
                -- them for clarity and efficiency
                if Row /= Ipass then
                   Factor := A(Row, Ipass);
                   for Col in Internal_Index_Type loop
                      if Row /= Ipass then
                         Inverse(Row, Col) := Inverse(Row, Col) - Factor * A(Ipass, Col);
                         A(Row, Col) := A(Row, Col) - Factor * A(Ipass, Col);
                      end if;
                   end loop;
                end if;
             end loop;

          end loop; -- iPass

       end Internal_Matrix_Inversion;
       pragma Inline (Internal_Matrix_Inversion);
    begin
       if A'Length(1)       = A'Length(2)       and
          Inverse'Length(1) = Inverse'Length(2) and
          A'Length(1)       = Inverse'Length(2)
       then
          Internal_Matrix_Inversion(A, Inverse);
       else
          raise Constraint_Error;
       end if;
    end Matrix_Inversion;

    -- #define ns 30
    -- #define ms 12
    Rs : constant := 6;
    -- A is an n by n matrix, A is the state transistion matrix
    -- B is an n by m matrix, B is the input matrix
    -- C is an r by n matrix, C is the measurement matrix
    -- xhat is an n by 1 vector, xhat is the current state of the system
    -- y is an r by 1 vector, y is the measured output of the system
    -- u is an m by 1 vector, u is the known input to the system
    -- Sz is an r by r matrix, Sz is the measurement nise covariance
    -- Sw is an n by n matrix, Sw is the process noise covariance
    -- P is an n by n matrix, P is the estimator error covariance
    -- float A[ns][ns];
    -- float B[ns][ms];
    -- float C[rs][ns];
    -- float xhat[ns][1];
    -- float y[rs][1];
    -- float u[ms][1];
    -- float Sz[rs][rs];
    -- float Sw[ns][ns];
    -- float P[ns][ns];
    -- float AP[ns][ns]; // This is the matrix A*P
    -- float CT[ns][rs]; // This is the matrix CT
    -- float APCT[ns][rs]; // This is the matrix A*P*CT
    -- float CP[rs][ns]; // This is the matrix C*P
    -- float CPCT[rs][rs]; // This is the matrix C*P*CT
    Cpctsz : Matrices (1 .. Rs, 1 .. Rs); -- This is the matrix C*P*CT+Sz
    --  CPCTSzInv: Matrices(1 .. Rs, 1 .. Rs); -- (C*P*CT+Sz)**-1
    Cpctszinv : Matrices (0 .. Rs - 1, 0 .. Rs - 1); -- (C*P*CT+Sz)**-1
    -- float K[ns][rs]; // This is the Kalman gain.
    -- float Cxhat[rs][1]; // This is the vector C*xhat
    -- float yCxhat[rs][1]; // This is the vector y-C*xhat
    -- float KyCxhat[ns][1]; // This is the vector K*(y-C*xhat)
    -- float Axhat[ns][1]; // This is the vector A*xhat
    -- float Bu[ns][1]; // This is the vector B*u
    -- float AxhatBu[ns][1]; // This is the vector A*xhat+B*u
    -- float AT[ns][ns]; // This is the matrix AT
    -- float APAT[ns][ns]; // This is the matrix A*P*AT
    -- float APATSw[ns][ns]; // This is the matrix A*P*AT+Sw
    -- float CPAT[rs][ns]; // This is the matrix C*P*AT
    -- float SzInv[rs][rs]; // This is the matrix Sz-1
    -- float APCTSzInv[ns][rs]; // This is the matrix A*P*CT*Sz-1
    -- float APCTSzInvCPAT[ns][ns]; // This is the matrix A*P*CT*Sz-1*C*P*AT
begin
    Matrix_Inversion(Cpctsz, Cpctszinv);
end Inv2;


Michael