C to Ada conversion

C to Ada conversion

[Mark]

Interested in doing some benchmarking in Ada.  Trouble is the I need
to convert the C code below to Ada and my Ada is well...
I'd greatly appreaciate it if someone could 'convert' the function
referenced below.  I'm also curious on how Ada allocates and
deallocates objects

Thanks 


--------------------------------------------------------------------

#define ns 30
#define ms 12
#define rs 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
float CPCTSz[rs][rs]; // This is the matrix C*P*CT+Sz
float CPCTSzInv[rs][rs]; // This is the matrix (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


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.
  int i, j, iPass, imx, icol, irow;
  float det, temp, pivot, factor;
  float* ac = (float*)calloc(n*n, sizeof(float));
  det = 1;
  for (i = 0; i < n; i++)
	{
	for (j = 0; j < n; j++)
		{
		AInverse[n*i+j] = 0;
		ac[n*i+j] = A[n*i+j];
		}
	AInverse[n*i+i] = 1;
	}
  // The current pivot row is iPass.  

  // For each pass, first find the maximum element in the pivot
column.
  for (iPass = 0; iPass < n; iPass++)
	{
	imx = iPass;
	for (irow = iPass; irow < n; irow++)
		{
		if (fabs(A[n*irow+iPass]) > fabs(A[n*imx+iPass])) imx = irow;
		}
	// Interchange the elements of row iPass and row imx in both A and
AInverse.
	if (imx != iPass)
		{
		for (icol = 0; icol < n; icol++)
			{
			temp = AInverse[n*iPass+icol];
			AInverse[n*iPass+icol] = AInverse[n*imx+icol];
			AInverse[n*imx+icol] = temp;
			if (icol >= iPass)
				{
				temp = A[n*iPass+icol];
				A[n*iPass+icol] = A[n*imx+icol];
				A[n*imx+icol] = temp;
				}
			}
		}
	// The current pivot is now A[iPass][iPass].
	// The determinant is the product of the pivot elements.
	pivot = A[n*iPass+iPass];
	det = det * pivot;
	if (det == 0) 
		{
		free(ac);
		return 0;
		}
	  for (icol = 0; icol < n; icol++)
		{
		// Normalize the pivot row by dividing by the pivot element.
		AInverse[n*iPass+icol] = AInverse[n*iPass+icol] / pivot;
		if (icol >= iPass) A[n*iPass+icol] = A[n*iPass+icol] / pivot;
		}
	  for (irow = 0; irow < n; irow++)
		// 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.
		{
		if (irow != iPass) factor = A[n*irow+iPass];
		for (icol = 0; icol < n; icol++)
			{
			if (irow != iPass)
				{
				AInverse[n*irow+icol] -= factor * AInverse[n*iPass+icol];
				A[n*irow+icol] -= factor * A[n*iPass+icol];
				}
			}
		}
	}
	free(ac);
	return 1;
}

int main (void)
{

  int n = ns;
  int m = ms;
  int r = rs;

  MatrixInversion((float*)CPCTSz, r, (float*)CPCTSzInv);
  return 0;
}

Re: C to Ada conversion

[tmoran]

> Interested in doing some benchmarking in Ada.  ... convert
  See below
> I'm also curious on how Ada allocates and deallocates objects
  There is an alloc(new)/free mechanism but usually it's more efficient,
as well as easier and less error prone, to just let the compiler allocate
things stack-wise.  See matrix "Ac" below.  It's allocated automatically,
and deallocated automatically.  (Actually, since it's never used here,
one hopes a good compiler doesn't bother generating any code for it.)

procedure Inv is

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

  procedure Matrix_Inversion(A       : in out Matrices;
                             Inverse :    out 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 Inverse'range (1) loop
      Inverse(I, I) := 1.0;
    end loop;

    for iPass in A'range (1) loop
      -- For each pass, first find maximum element in the pivot column.
      declare
        Max_Row : Integer := iPass;
      begin
        for Row in iPass .. A'Last(1) 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 A'range (2) 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 A'range (2) 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 A'range (1) 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 A'range (2) 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 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
  -- 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 Inv;

Re: C to Ada conversion

[Jeffrey Carter]

tmoran@acm.org wrote:
>>Interested in doing some benchmarking in Ada.  ... convert
> 
>   See below

Since the original is from UCF, I suspect you've just done his homework 
for him.

-- 
Jeff Carter
"I'm a lumberjack and I'm OK."
Monty Python's Flying Circus

Re: C to Ada conversion

[tmoran]

> Since the original is from UCF, I suspect you've just done his homework
> for him.
  I hope that's not the case, or if it is, that his teacher has high
expectations for the rest of the semester. ;)

Re: C to Ada conversion

[Mark]

tmoran@acm.org wrote in message news:<giLU9.550780$GR5.330948@rwcrnsc51.ops.asp.att.net>...
> > Since the original is from UCF, I suspect you've just done his homework
> > for him.
>   I hope that's not the case, or if it is, that his teacher has high
> expectations for the rest of the semester. ;)


Homework problem.  Oh no!!!  The harsh reality of it is 'Universities'
promote C and it's cousin ++.  That said, it's highly unlikely a
professor would ask a student 'today' to convert C code to Ada.
The UCF account is my 'better' email .. so to speak.

Re: C to Ada conversion

[John R. Strohm]

"Jeffrey Carter" <jrcarter@acm.org> wrote in message
news:3E2367A1.5050107@acm.org...
> tmoran@acm.org wrote:
> >>Interested in doing some benchmarking in Ada.  ... convert
> >
> >   See below
>
> Since the original is from UCF, I suspect you've just done his homework
> for him.

Unlikely.  Most professors don't assign Kalman filters on homework
assignments.

Re: C to Ada conversion

[Martin Dowie]

<tmoran@acm.org> wrote in message news:ELHU9.51517$%n.18940@sccrnsc02...
> > Interested in doing some benchmarking in Ada.  ... convert
>   See below
[snip]

Did you perform any benchmarking between the two versions?

Re: C to Ada conversion

[tmoran]

> Did you perform any benchmarking between the two versions?
  Done right, that's non trivial.

Re: C to Ada conversion

[Martin Dowie]

<tmoran@acm.org> wrote in message news:zaQU9.674596$WL3.714902@rwcrnsc54...
> > Did you perform any benchmarking between the two versions?
>   Done right, that's non trivial.

Absolutely - especially in a Win32-type environment with other
processes coming in and out.

I was also wondering what differences might show up between
compilers - both C and Ada.

Re: C to Ada conversion

[Mark]

Martin I've done some benchmarking in C and quite frankly the ONLY
reason I'm attempting the Ada benchmark is to get a feel for how the
Ada code will perform on an existing system and get a feel for some of
the 'fuss' I hear about how Ada is verbose, etc....

"Martin Dowie" <martin.dowie@no.spam.btopenworld.com> wrote in message news:<b00dnn$pj8$1@venus.btinternet.com>...
> <tmoran@acm.org> wrote in message news:ELHU9.51517$%n.18940@sccrnsc02...
> > > Interested in doing some benchmarking in Ada.  ... convert
> >   See below
> [snip]
> 
> Did you perform any benchmarking between the two versions?

Re: C to Ada conversion

[Dr. Michael Paus]

tmoran@acm.org wrote:
>>Interested in doing some benchmarking in Ada.  ... convert
> 
>   See below
> 
>>I'm also curious on how Ada allocates and deallocates objects
> 
>   There is an alloc(new)/free mechanism but usually it's more efficient,
> as well as easier and less error prone, to just let the compiler allocate
> things stack-wise.  See matrix "Ac" below.  It's allocated automatically,
> and deallocated automatically.  (Actually, since it's never used here,
> one hopes a good compiler doesn't bother generating any code for it.)
[...]

I hope nobody will be using this Ada code for any serious purpose.
As a quick and dirty conversion for some benchmark it is certainly
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
start with 1 (or 0 as in C). A general purpose matrix package should
be able to cope with the situation that the index range could be
anything as long as the involved matrices are compatibel in size.
Otherwise it would be difficult to use slices of arrays for example.

If someone uses the example here with a definition like this

   CPCTSzInv: Matrices(0 .. Rs - 1, 0 .. Rs - 1); -- (C*P*CT+Sz)**-1

with the other matrix unchanged the compiler won't complain but the
result will be a constraint_error at run-time. In C you simply don't
have this specific Ada problem because an array by definition has
an index range starting with 0 and there is no good reason why the
inversion should not be possible in this case. In practice you often
have the need to mix arrays which have compatibel sizes but not the
same index range.

Just my 2 (Euro-)cents

Michael

Re: C to Ada conversion

[tmoran]

>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
>start with 1 (or 0 as in C). A general purpose matrix package should
>be able to cope with the situation that the index range could be
>anything as long as the involved matrices are compatibel in size.
  That's why the comments said:
>>  -- 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.
  For quick transliteration that didn't seem important.  For a benchmark
to compare with C, the extra Ada generality seemed like a bad idea.  For
a general purpose matrix package, I agree, the considerations are
different.

Re: C to Ada conversion

[tmoran]

> >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;

Re: C to Ada conversion

[Dr. Michael Paus]

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

Re: C to Ada conversion

[Vinzent Hoefler]

Dr. Michael Paus wrote:

>    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.

Turning off checks does not mean that Constraint_Error will not be 
raised (especially if you consider the divide-by-zero case).


Vinzent.

-- 
A general leading the State Department resembles  a dragon commanding
ducks.
                -- New York Times, Jan. 20, 1981

Re: C to Ada conversion

[tmoran]

> >   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)
  Oops.  I saw the #define on rs before the array declarations and thought
"Aha, known at compile time", but of course you are right that the
Matrix_Inversion routine doesn't depend on rs and will take any size
array.  So I withdraw my suggestion.

> This is NOT pointless because a piece of code should also
> work correctly if someone turns off the automatic checks in Ada.
  I disagree.  Turning off checks is a substantial change to the
assumptions of the programmer.  If he coded to not depend on the
assumption that checks are on, great.  But I don't think one should
assume the programmer didn't assume checks unless he says so.  If
the code should be able to run with checks off, there's no reason
to ever have them on, and the source code should include the
pragma to turn them off.

  Converting index ranges like this by passing to a subroutine is
very nice.  I hadn't seen that done before.

Re: C to Ada conversion

[Mark]

Gents, Could someone point me in the right direction on converting a
few other routines.  Similar to that done by tmorgan.  Dont think i
did things right.  These are much simpler.

Here's my definitons

  Ns : constant Integer := 30;
  Ms : constant Integer := 12;
  Rs : constant Integer := 6;

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

  -- 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

  A         : Matrices(1 .. Ns, 1 .. Ns);
  B         : Matrices(1 .. Ns, 1 .. Ms);
  C         : Matrices(1 .. Ns, 1 .. Rs);
  xhat      : Matrices(1 .. Ns, 1 .. 1);
  Y         : Matrices(1 .. Rs, 1 .. 1);
  U         : Matrices(1 .. Ms, 1 .. 1);
  Sz        : Matrices(1 .. Rs, 1 .. Rs);
  Sw        : Matrices(1 .. Ns, 1 .. Ns);
  P         : Matrices(1 .. Ns, 1 .. Ns);
  AP        : Matrices(1 .. Ns, 1 .. Ns);  -- This is the matrix A*P
  CT        : Matrices(1 .. Ns, 1 .. Rs);  -- This is the matrix CT
  APCT      : Matrices(1 .. Ns, 1 .. Rs);  -- This is the matrix
A*P*CT
  CP        : Matrices(1 .. Rs, 1 .. Ns);  -- This is the matrix C*P
  CPCT      : Matrices(1 .. Rs, 1 .. 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
  K         : Matrices(1 .. Ns, 1 .. Rs);  -- This is the Kalman gain.
  Cxhat     : Matrices(1 .. Rs, 1 .. 1);   -- This is the vector
C*xhat
  yCxhat    : Matrices(1 .. Rs, 1 .. 1);   -- This is the vector
y-C*xhat
  KyCxhat   : Matrices(1 .. Ns, 1 .. 1);   -- This is the vector
K*(y-C*xhat)
  Axhat     : Matrices(1 .. Ns, 1 .. 1);   -- This is the vector
A*xhat
  Bu        : Matrices(1 .. Ns, 1 .. 1);   -- This is the vector B*u
  AxhatBu   : Matrices(1 .. Ns, 1 .. 1);   -- This is the vector
A*xhat+B*u
  AT        : Matrices(1 .. Ns, 1 .. Ns);  -- This is the matrix AT
  APAT      : Matrices(1 .. Ns, 1 .. Ns);  -- This is the matrix
A*P*AT
  APATSw    : Matrices(1 .. Ns, 1 .. Ns);  -- This is the matrix
A*P*AT+Sw
  CPAT      : Matrices(1 .. Rs, 1 .. Ns);  -- This is the matrix
C*P*AT
  SzInv     : Matrices(1 .. Rs, 1 .. Rs);  -- This is the matrix Sz-1
  APCTSzInv : Matrices(1 .. Ns, 1 .. Rs);  -- This is the matrix
A*P*CT*Sz-1
  APCTSzInvCPAT  : Matrices(1 .. Ns, 1 .. Ns);   -- This is the matrix
A*P*CT*Sz-1*C*P*AT


Thanks in advance. 


void MatrixMultiply(float* A, float* B, int m, int p, int n, float* C)
{       

  /* 
  A = input matrix (m x p)
  B = input matrix (p x n)
  m = number of rows in A
  p = number of columns in A = number of columns in B       
  n = number of columns in B
  C = output matrix = A*B (m x n)
  */ 
  int i, j, k;
  for (i=0;i<m;i++)
    for(j=0;j<n;j++)
	{
	  C[n*i+j]=0;
	  for (k=0;k<p;k++)
	    C[n*i+j]= C[n*i+j]+A[p*i+k]*B[n*k+j];
	}
}

-- following tmorgans convesion
procedure MatrixMultiply(A: in out Matrices;
                         
begin

end MatrixMultiply;

-----
void MatrixAddition(float* A, float* B, int m, int n, float* C)
{
  /* 
  A = input matrix (m x n)
  B = input matrix (m x n)
  m = number of rows in A = number of rows in B
  n = number of columns in A = number of columns in B
  C = output matrix = A+B (m x n)  
  */ 
  int i, j;
  for (i=0;i<m;i++)
    for(j=0;j<n;j++)
	  C[n*i+j]=A[n*i+j]+B[n*i+j];
}
void MatrixSubtraction(float* A, float* B, int m, int n, float* C)
{
  /* 
  A = input matrix (m x n)
  B = input matrix (m x n)
  m = number of rows in A = number of rows in B
  n = number of columns in A = number of columns in B
  C = output matrix = A-B (m x n)  
  */
  int i, j;
  for (i=0;i<m;i++)
    for(j=0;j<n;j++)
	  C[n*i+j]=A[n*i+j]-B[n*i+j];
}                    

void MatrixTranspose(float* A, int m, int n, float* C)
{    

     /* 
     A = input matrix (m x n)
     m = number of rows in A 
     n = number of columns in A 
     C = output matrix = the transpose of A (n x m)
     */
  int i, j;
  for (i=0;i<m;i++)
    for(j=0;j<n;j++)
      C[m*j+i]=A[n*i+j];     
      
}


tmoran@acm.org wrote in message news:<rH%U9.74545$3v.13364@sccrnsc01>...
> > >   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)
>   Oops.  I saw the #define on rs before the array declarations and thought
> "Aha, known at compile time", but of course you are right that the
> Matrix_Inversion routine doesn't depend on rs and will take any size
> array.  So I withdraw my suggestion.
> 
> > This is NOT pointless because a piece of code should also
> > work correctly if someone turns off the automatic checks in Ada.
>   I disagree.  Turning off checks is a substantial change to the
> assumptions of the programmer.  If he coded to not depend on the
> assumption that checks are on, great.  But I don't think one should
> assume the programmer didn't assume checks unless he says so.  If
> the code should be able to run with checks off, there's no reason
> to ever have them on, and the source code should include the
> pragma to turn them off.
> 
>   Converting index ranges like this by passing to a subroutine is
> very nice.  I hadn't seen that done before.

Re: C to Ada conversion

[tmoran]

Do you want examples of conversion, in which case doing both add and
subtract is surely redundant, or do you want a library of matrix
operations?  I've gotten the latter from the Numerical Analysis
Group and the Walnut Creek or ASE2 CDROMs, and I'm sure other places.

Re: C to Ada conversion

[John R. Strohm]

Y'know, if you asked around politely, you could probably find a
completely-written Kalman filter in Ada somewhere.  I'd try looking in the
Aerospace Engineering department: some of the professors probably know
people at the big defense companies.

"Mark" <ma740988@pegasus.cc.ucf.edu> wrote in message
news:a5ae824.0301161317.3e084a77@posting.google.com...
> Gents, Could someone point me in the right direction on converting a
> few other routines.  Similar to that done by tmorgan.  Dont think i
> did things right.  These are much simpler.
>
> Here's my definitons
>
>   Ns : constant Integer := 30;
>   Ms : constant Integer := 12;
>   Rs : constant Integer := 6;
>
>   type Matrices is array (Integer range <>, Integer range <>) of
> Float;
>
>   -- 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
>
>   A         : Matrices(1 .. Ns, 1 .. Ns);
>   B         : Matrices(1 .. Ns, 1 .. Ms);
>   C         : Matrices(1 .. Ns, 1 .. Rs);
>   xhat      : Matrices(1 .. Ns, 1 .. 1);
>   Y         : Matrices(1 .. Rs, 1 .. 1);
>   U         : Matrices(1 .. Ms, 1 .. 1);
>   Sz        : Matrices(1 .. Rs, 1 .. Rs);
>   Sw        : Matrices(1 .. Ns, 1 .. Ns);
>   P         : Matrices(1 .. Ns, 1 .. Ns);
>   AP        : Matrices(1 .. Ns, 1 .. Ns);  -- This is the matrix A*P
>   CT        : Matrices(1 .. Ns, 1 .. Rs);  -- This is the matrix CT
>   APCT      : Matrices(1 .. Ns, 1 .. Rs);  -- This is the matrix
> A*P*CT
>   CP        : Matrices(1 .. Rs, 1 .. Ns);  -- This is the matrix C*P
>   CPCT      : Matrices(1 .. Rs, 1 .. 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
>   K         : Matrices(1 .. Ns, 1 .. Rs);  -- This is the Kalman gain.
>   Cxhat     : Matrices(1 .. Rs, 1 .. 1);   -- This is the vector
> C*xhat
>   yCxhat    : Matrices(1 .. Rs, 1 .. 1);   -- This is the vector
> y-C*xhat
>   KyCxhat   : Matrices(1 .. Ns, 1 .. 1);   -- This is the vector
> K*(y-C*xhat)
>   Axhat     : Matrices(1 .. Ns, 1 .. 1);   -- This is the vector
> A*xhat
>   Bu        : Matrices(1 .. Ns, 1 .. 1);   -- This is the vector B*u
>   AxhatBu   : Matrices(1 .. Ns, 1 .. 1);   -- This is the vector
> A*xhat+B*u
>   AT        : Matrices(1 .. Ns, 1 .. Ns);  -- This is the matrix AT
>   APAT      : Matrices(1 .. Ns, 1 .. Ns);  -- This is the matrix
> A*P*AT
>   APATSw    : Matrices(1 .. Ns, 1 .. Ns);  -- This is the matrix
> A*P*AT+Sw
>   CPAT      : Matrices(1 .. Rs, 1 .. Ns);  -- This is the matrix
> C*P*AT
>   SzInv     : Matrices(1 .. Rs, 1 .. Rs);  -- This is the matrix Sz-1
>   APCTSzInv : Matrices(1 .. Ns, 1 .. Rs);  -- This is the matrix
> A*P*CT*Sz-1
>   APCTSzInvCPAT  : Matrices(1 .. Ns, 1 .. Ns);   -- This is the matrix
> A*P*CT*Sz-1*C*P*AT
>
>
> Thanks in advance.
>
>
> void MatrixMultiply(float* A, float* B, int m, int p, int n, float* C)
> {
>
>   /*
>   A = input matrix (m x p)
>   B = input matrix (p x n)
>   m = number of rows in A
>   p = number of columns in A = number of columns in B
>   n = number of columns in B
>   C = output matrix = A*B (m x n)
>   */
>   int i, j, k;
>   for (i=0;i<m;i++)
>     for(j=0;j<n;j++)
> {
>   C[n*i+j]=0;
>   for (k=0;k<p;k++)
>     C[n*i+j]= C[n*i+j]+A[p*i+k]*B[n*k+j];
> }
> }
>
> -- following tmorgans convesion
> procedure MatrixMultiply(A: in out Matrices;
>
> begin
>
> end MatrixMultiply;
>
> -----
> void MatrixAddition(float* A, float* B, int m, int n, float* C)
> {
>   /*
>   A = input matrix (m x n)
>   B = input matrix (m x n)
>   m = number of rows in A = number of rows in B
>   n = number of columns in A = number of columns in B
>   C = output matrix = A+B (m x n)
>   */
>   int i, j;
>   for (i=0;i<m;i++)
>     for(j=0;j<n;j++)
>   C[n*i+j]=A[n*i+j]+B[n*i+j];
> }
> void MatrixSubtraction(float* A, float* B, int m, int n, float* C)
> {
>   /*
>   A = input matrix (m x n)
>   B = input matrix (m x n)
>   m = number of rows in A = number of rows in B
>   n = number of columns in A = number of columns in B
>   C = output matrix = A-B (m x n)
>   */
>   int i, j;
>   for (i=0;i<m;i++)
>     for(j=0;j<n;j++)
>   C[n*i+j]=A[n*i+j]-B[n*i+j];
> }
>
> void MatrixTranspose(float* A, int m, int n, float* C)
> {
>
>      /*
>      A = input matrix (m x n)
>      m = number of rows in A
>      n = number of columns in A
>      C = output matrix = the transpose of A (n x m)
>      */
>   int i, j;
>   for (i=0;i<m;i++)
>     for(j=0;j<n;j++)
>       C[m*j+i]=A[n*i+j];
>
> }
>
>
> tmoran@acm.org wrote in message news:<rH%U9.74545$3v.13364@sccrnsc01>...
> > > >   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)
> >   Oops.  I saw the #define on rs before the array declarations and
thought
> > "Aha, known at compile time", but of course you are right that the
> > Matrix_Inversion routine doesn't depend on rs and will take any size
> > array.  So I withdraw my suggestion.
> >
> > > This is NOT pointless because a piece of code should also
> > > work correctly if someone turns off the automatic checks in Ada.
> >   I disagree.  Turning off checks is a substantial change to the
> > assumptions of the programmer.  If he coded to not depend on the
> > assumption that checks are on, great.  But I don't think one should
> > assume the programmer didn't assume checks unless he says so.  If
> > the code should be able to run with checks off, there's no reason
> > to ever have them on, and the source code should include the
> > pragma to turn them off.
> >
> >   Converting index ranges like this by passing to a subroutine is
> > very nice.  I hadn't seen that done before.