Re: C to Ada conversion

["John R. Strohm" <strohm@airmail.net>]

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.