-- TEST_LUD_PACKAGE.ADA LU decomposition -- Decompose a matrix A into a Lower triangular matrix -- times an Upper triangluar matrix A = L * U -- where the diagonal of L is all ones. -- Or more generally P * A = L * U where P is a permutation -- matrix (a permutation of the rows of the Identity matrix) -- represented as an integer vector of subscripts PV=2,3,1,4 -- means a 1.0 in P at 1,2 2,3 3,1 4,4 (P**t = P**-1) -- A = P**T * L * U -- -- Array_Index_Error is raised if matrix not square and also if -- A, L and U row and column ranges are not the same. -- -- -- programmed 3/5/94 by Jon S. Squire with LUD_PACKAGE; use LUD_PACKAGE; with REAL_ARRAYS; use REAL_ARRAYS; with REAL_ARRAYS_IO; with INTEGER_ARRAYS; use INTEGER_ARRAYS; with INTEGER_ARRAYS_IO; with TEXT_IO; use TEXT_IO; procedure TEST_LUD_PACKAGE is subtype real is float; package real_arr_io is new real_arrays_io(real, real_vector, real_matrix); use real_arr_io; package int_arr_io is new integer_arrays_io (integer, integer_vector, integer_matrix); use int_arr_io; package int_io is new integer_io(integer); use int_io; n : integer; m : integer; Input : text_io.File_TYpe; Input_File_Name : string(1..260); -- "-" means read from Standard Input Input_file_name_length : Natural; Title : string(1..260); Title_length : natural; debug : boolean := false; -- turned on by $ on front of Input_file_name finished : exception; array_index_error : exception renames Real_Arrays.Array_Index_Error; function norm1(A : Real_Matrix) return real is value : real := 0.0; begin for i in A'range(1) loop for j in A'range(2) loop if abs A(i,j) > value then value := abs A(i,j); end if; end loop; end loop; return value; end norm1; function norm1(A : Real_vector) return real is value : real := 0.0; begin for i in A'range loop if abs A(i) > value then value := abs A(i); end if; end loop; return value; end norm1; function MAKE_PERM( P : INTEGER_VECTOR ) return REAL_MATRIX is PP : REAL_MATRIX(P'RANGE,P'RANGE) := (others=>(others=>0.0)); begin for I in P'RANGE loop PP(I,P(I)) := 1.0; end loop; return PP; end MAKE_PERM; begin put_line("test_lud_package, enter data file name:"); get_line(input_file_name, input_file_name_length); if input_file_name_length = 0 then put_line("no file name provided, ending."); raise finished; end if; put_line(input_file_name(1..input_file_name_length)); if input_file_name(1) = '$' then debug := true; input_file_name(1..input_file_name_length-1) := input_file_name(2..input_file_name_length); input_file_name_length := input_file_name_length-1; end if; -- test for '-' open(input, in_file, input_file_name(1..input_file_name_length)); loop -- read data sets from input file, stop on end of file or "0" -- get n and m int_io.get(Input,n); if n = 0 then raise finished; end if; int_io.get(Input,m); get_line(Input, Title, Title_length); put_line(" n=" & integer'image(n) & " m=" & integer'image(m) & " " & Title(1..Title_length)); declare A : real_matrix(1..m,1..n); AA : real_matrix(1..m,1..n); AI : real_matrix(1..n,1..m); L : real_matrix(1..m,1..n); IDENT : real_matrix(1..n,1..n); U : real_matrix(1..n,1..n); P : integer_vector(1..n); PP : real_matrix(1..n,1..n); X : real_vector(1..n); Y : real_vector(1..n); Z : real_vector(1..n); sum : real; begin get(Input,A); -- read in matrix for i in Y'range loop Y(i) := real(i); end loop; if debug then put_line("A= "); put(A); end if; lud(A, L, U, P); if debug then put_line("L= "); put(L); put_line("U= "); put(U); put_line("P= "); put(P); new_line; end if; PP := MAKE_PERM(P); AA := PP * A - L * U; put(norm1(AA)); put_line("= norm1 of PP * A - L * U"); AA := TRANSPOSE(PP) * L * U ; put(norm1(AA-A)); put_line("= norm1 of A - (transpose(PP) * L * U)"); sum := 0.0; for I in U'first(1)+1..U'last(1) loop for J in U'first(2)..I-1-U'first(1)+U'first(2) loop sum := sum + abs U(I,J); end loop; end loop; put(sum); put_line("= norm1 of U below diagonal"); sum := 0.0; for I in L'first(1)..L'last(1) loop sum := sum + abs(1.0-L(I,I-L'first(1)+L'first(2))); for J in I+1-L'first(1)+L'first(2)..L'last(2) loop sum := sum + abs L(I,J); end loop; end loop; put(sum); put_line("= norm1 of L-I on and above diagonal"); if debug then put_line("AA=transpose(PP)*L*U"); put(AA); new_line; end if; X := lud_solve(L, U, P, Y); Z := A * X; if debug then put_line("Y= "); put(Y); put_line("X= "); put(X); put_line("Z= "); put(Z); put_line("Z should equal Y, based on A"); new_line; end if; put(norm1(Z-Y)); put_line("= norm1 of origonal matrix * X - Y, solving equations"); AI := lud_inverse(L, U, P); if debug then put_line("lud inverse AI"); put(AI); new_line; put_line("A * AI should be identity"); IDENT := A * AI; put(IDENT); new_line; put_line("AI * A "); IDENT := AI * A; put(IDENT); new_line; end if; put(norm1(A*AI-Identity_Matrix(n))); put_line("= norm1 of origonal matrix * inverse - identity"); put(norm1(AI*A-Identity_Matrix(n))); put_line("= norm1 of inverse * origonal matrix - identity"); put(norm1(AA*AI-Identity_Matrix(n))); put_line("= norm1 of AA * AI - I"); exception when Array_Index_error => put_line("Array_Index_Error from package"); when Singular_error => put_line("singular matrix"); new_line; end; -- tests on each data set end loop; -- all data sets finished put_line("test_lud_package completed"); exception when Text_io.End_Error => put_line("test_lud_package completed"); when DATA_ERROR | LAYOUT_ERROR => put_line("bad input data, terminating at:"); put_line(Title(1..Title_length)); when finished => put_line("test_lud_package ends"); when STATUS_ERROR | MODE_ERROR | NAME_ERROR | USE_ERROR | DEVICE_ERROR => put_line("Could not open: " & input_file_name(1..input_file_name_length)); -- when others => -- put_line("FAILED premature exit of test by unhandled exception"); end TEST_LUD_PACKAGE;