From mboxrd@z Thu Jan 1 00:00:00 1970 X-Spam-Checker-Version: SpamAssassin 3.4.4 (2020-01-24) on polar.synack.me X-Spam-Level: X-Spam-Status: No, score=-1.9 required=5.0 tests=BAYES_00 autolearn=ham autolearn_force=no version=3.4.4 X-Google-Language: ENGLISH,ASCII-7-bit X-Google-Thread: 103376,83d2d63f98e99c58 X-Google-Attributes: gid103376,public X-Google-ArrivalTime: 2003-01-14 11:48:29 PST Path: archiver1.google.com!news1.google.com!newsfeed.stanford.edu!bloom-beacon.mit.edu!npeer.de.kpn-eurorings.net!rz.uni-karlsruhe.de!schlund.de!news.online.de!not-for-mail From: "Dr. Michael Paus" Newsgroups: comp.lang.ada Subject: Re: C to Ada conversion Date: Tue, 14 Jan 2003 20:48:22 +0100 Organization: 1&1 Internet AG Message-ID: References: NNTP-Posting-Host: p5083059c.dip0.t-ipconnect.de Mime-Version: 1.0 Content-Type: text/plain; charset=us-ascii; format=flowed Content-Transfer-Encoding: 7bit X-Trace: news.online.de 1042573703 31539 80.131.5.156 (14 Jan 2003 19:48:23 GMT) X-Complaints-To: abuse@online.de NNTP-Posting-Date: 14 Jan 2003 19:48:23 GMT User-Agent: Mozilla/5.0 (Windows; U; Windows NT 5.1; en-US; rv:1.2.1) Gecko/20021130 X-Accept-Language: en-us, en In-Reply-To: Xref: archiver1.google.com comp.lang.ada:33011 Date: 2003-01-14T19:48:23+00:00 List-Id: 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