* C to Ada conversion @ 2003-01-13 21:22 Mark 2003-01-13 23:14 ` tmoran 0 siblings, 1 reply; 19+ messages in thread From: Mark @ 2003-01-13 21:22 UTC (permalink / raw) 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; } ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-13 21:22 C to Ada conversion Mark @ 2003-01-13 23:14 ` tmoran 2003-01-14 1:27 ` Jeffrey Carter ` (2 more replies) 0 siblings, 3 replies; 19+ messages in thread From: tmoran @ 2003-01-13 23:14 UTC (permalink / raw) > 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; ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-13 23:14 ` tmoran @ 2003-01-14 1:27 ` Jeffrey Carter 2003-01-14 3:16 ` tmoran 2003-01-14 19:35 ` John R. Strohm 2003-01-14 7:20 ` Martin Dowie 2003-01-14 17:18 ` Dr. Michael Paus 2 siblings, 2 replies; 19+ messages in thread From: Jeffrey Carter @ 2003-01-14 1:27 UTC (permalink / raw) 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 ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 1:27 ` Jeffrey Carter @ 2003-01-14 3:16 ` tmoran 2003-01-14 12:52 ` Mark 2003-01-14 19:35 ` John R. Strohm 1 sibling, 1 reply; 19+ messages in thread From: tmoran @ 2003-01-14 3:16 UTC (permalink / raw) > 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. ;) ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 3:16 ` tmoran @ 2003-01-14 12:52 ` Mark 0 siblings, 0 replies; 19+ messages in thread From: Mark @ 2003-01-14 12:52 UTC (permalink / raw) 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. ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 1:27 ` Jeffrey Carter 2003-01-14 3:16 ` tmoran @ 2003-01-14 19:35 ` John R. Strohm 1 sibling, 0 replies; 19+ messages in thread From: John R. Strohm @ 2003-01-14 19:35 UTC (permalink / raw) "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. ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-13 23:14 ` tmoran 2003-01-14 1:27 ` Jeffrey Carter @ 2003-01-14 7:20 ` Martin Dowie 2003-01-14 8:49 ` tmoran 2003-01-14 13:01 ` Mark 2003-01-14 17:18 ` Dr. Michael Paus 2 siblings, 2 replies; 19+ messages in thread From: Martin Dowie @ 2003-01-14 7:20 UTC (permalink / raw) <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? ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 7:20 ` Martin Dowie @ 2003-01-14 8:49 ` tmoran 2003-01-14 22:45 ` Martin Dowie 2003-01-14 13:01 ` Mark 1 sibling, 1 reply; 19+ messages in thread From: tmoran @ 2003-01-14 8:49 UTC (permalink / raw) > Did you perform any benchmarking between the two versions? Done right, that's non trivial. ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 8:49 ` tmoran @ 2003-01-14 22:45 ` Martin Dowie 0 siblings, 0 replies; 19+ messages in thread From: Martin Dowie @ 2003-01-14 22:45 UTC (permalink / raw) <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. ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 7:20 ` Martin Dowie 2003-01-14 8:49 ` tmoran @ 2003-01-14 13:01 ` Mark 1 sibling, 0 replies; 19+ messages in thread From: Mark @ 2003-01-14 13:01 UTC (permalink / raw) 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? ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-13 23:14 ` tmoran 2003-01-14 1:27 ` Jeffrey Carter 2003-01-14 7:20 ` Martin Dowie @ 2003-01-14 17:18 ` Dr. Michael Paus 2003-01-14 18:10 ` tmoran 2 siblings, 1 reply; 19+ messages in thread From: Dr. Michael Paus @ 2003-01-14 17:18 UTC (permalink / raw) 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 ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 17:18 ` Dr. Michael Paus @ 2003-01-14 18:10 ` tmoran 2003-01-14 18:24 ` tmoran 0 siblings, 1 reply; 19+ messages in thread From: tmoran @ 2003-01-14 18:10 UTC (permalink / raw) >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. ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 18:10 ` tmoran @ 2003-01-14 18:24 ` tmoran 2003-01-14 19:48 ` Dr. Michael Paus 0 siblings, 1 reply; 19+ messages in thread From: tmoran @ 2003-01-14 18:24 UTC (permalink / raw) > >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; ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 18:24 ` tmoran @ 2003-01-14 19:48 ` Dr. Michael Paus 2003-01-14 20:12 ` Vinzent Hoefler 2003-01-14 21:55 ` tmoran 0 siblings, 2 replies; 19+ messages in thread From: Dr. Michael Paus @ 2003-01-14 19:48 UTC (permalink / raw) 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 ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 19:48 ` Dr. Michael Paus @ 2003-01-14 20:12 ` Vinzent Hoefler 2003-01-14 21:55 ` tmoran 1 sibling, 0 replies; 19+ messages in thread From: Vinzent Hoefler @ 2003-01-14 20:12 UTC (permalink / raw) 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 ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 19:48 ` Dr. Michael Paus 2003-01-14 20:12 ` Vinzent Hoefler @ 2003-01-14 21:55 ` tmoran 2003-01-16 21:17 ` Mark 1 sibling, 1 reply; 19+ messages in thread From: tmoran @ 2003-01-14 21:55 UTC (permalink / raw) > > 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. ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-14 21:55 ` tmoran @ 2003-01-16 21:17 ` Mark 2003-01-17 0:10 ` tmoran 2003-01-17 1:09 ` John R. Strohm 0 siblings, 2 replies; 19+ messages in thread From: Mark @ 2003-01-16 21:17 UTC (permalink / raw) 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. ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-16 21:17 ` Mark @ 2003-01-17 0:10 ` tmoran 2003-01-17 1:09 ` John R. Strohm 1 sibling, 0 replies; 19+ messages in thread From: tmoran @ 2003-01-17 0:10 UTC (permalink / raw) 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. ^ permalink raw reply [flat|nested] 19+ messages in thread
* Re: C to Ada conversion 2003-01-16 21:17 ` Mark 2003-01-17 0:10 ` tmoran @ 2003-01-17 1:09 ` John R. Strohm 1 sibling, 0 replies; 19+ messages in thread From: John R. Strohm @ 2003-01-17 1:09 UTC (permalink / raw) 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. ^ permalink raw reply [flat|nested] 19+ messages in thread
end of thread, other threads:[~2003-01-17 1:09 UTC | newest] Thread overview: 19+ messages (download: mbox.gz / follow: Atom feed) -- links below jump to the message on this page -- 2003-01-13 21:22 C to Ada conversion Mark 2003-01-13 23:14 ` tmoran 2003-01-14 1:27 ` Jeffrey Carter 2003-01-14 3:16 ` tmoran 2003-01-14 12:52 ` Mark 2003-01-14 19:35 ` John R. Strohm 2003-01-14 7:20 ` Martin Dowie 2003-01-14 8:49 ` tmoran 2003-01-14 22:45 ` Martin Dowie 2003-01-14 13:01 ` Mark 2003-01-14 17:18 ` Dr. Michael Paus 2003-01-14 18:10 ` tmoran 2003-01-14 18:24 ` tmoran 2003-01-14 19:48 ` Dr. Michael Paus 2003-01-14 20:12 ` Vinzent Hoefler 2003-01-14 21:55 ` tmoran 2003-01-16 21:17 ` Mark 2003-01-17 0:10 ` tmoran 2003-01-17 1:09 ` John R. Strohm
This is a public inbox, see mirroring instructions for how to clone and mirror all data and code used for this inbox