-- TEST_QR_LINEAR.ADA QR decomposition -- Decompose a matrix A into an orthonormal matrix Q -- times an Upper triangluar matrix R such that -- A = Q * R -- Q * Q**T = I and R is an upper diagonal matrix -- A x = b or Q * R x = b becomes Q**T * Q * R x = Q**T b -- or R x = Q**T b solve by back substitution -- A**-1 = set of x vectors where set of b is unit vectors -- -- Array_Index_Error is raised if matrix not square and also if -- A, Q and R row and column ranges are not the same. -- -- programmed 9/9/96 by Jon S. Squire with GENERIC_REAL_LINEAR_EQUATIONS; with REAL_ARRAYS; use REAL_ARRAYS; with REAL_ARRAYS_IO; with TEXT_IO; use TEXT_IO; procedure TEST_QR_LINEAR is subtype real is float; package arr_io is new real_arrays_io(real, real_vector, real_matrix); use arr_io; package int_io is new integer_io(integer); use int_io; package real_linear_equations is new generic_real_linear_equations( real, real_vector, real_matrix); use real_linear_equations; Array_Index_Error : exception renames real_linear_equations.array_index_error; 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; 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; begin put_line("TEST_QR_LINEAR, 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 -- leading $ on file name for debug 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 eof in loop 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 get(Input,n); if n <= 0 then raise finished; end if; 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); AAA : real_matrix(1..m,1..n); Q : real_matrix(1..m,1..n); IDENT : real_matrix(1..n,1..n); R : real_matrix(1..n,1..n); RR : real_matrix(1..n,1..n); X : real_vector(1..n); Y : real_vector(1..n); YY : real_vector(1..n); Z : real_vector(1..n); sum : real; INFO : integer; 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; qr_decomposition(A, Q, R); if debug then put_line("Q= "); put(Q); put_line("R= "); put(R); new_line; end if; AA := Q * R; put(norm1(AA-A)); put_line("= norm1 of origonal matrix minus Q * R"); sum := 0.0; for I in R'first(1)+1..R'last(1) loop for J in R'first(2)..I-1-R'first(1)+R'first(2) loop sum := sum + abs R(I,J); end loop; end loop; put(sum); put_line("= norm1 of R below diagonal"); if debug then put_line("AA=Q*R"); put(AA); new_line; end if; X := qr_solve(Q, R, 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 := qr_inverse(Q, R); if debug then put_line("qr 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"); if debug then put_line("transpose(Q) * Q = I "); IDENT := transpose(Q) * Q; put(IDENT); new_line; put_line("Q * teanspose(Q) = I "); IDENT := Q * transpose(Q); put(IDENT); new_line; end if; put(norm1(Q*transpose(Q)-Identity_matrix(n))); put_line("= norm1 of Q * transpose(Q) - I"); put(norm1(transpose(Q)*Q-Identity_matrix(n))); put_line("= norm1 of transpose(Q) * Q - I"); exception when Array_Index_error => put_line("Array_Index_Error from package"); when Matrix_Data_Error => put_line("singular matrix"); new_line; end; -- tests on each data set end loop; -- all data sets finished put_line("TEST_QR_LINEAR completed"); exception when Text_io.End_Error => put_line("TEST_QR_LINEAR completed"); when DATA_ERROR | LAYOUT_ERROR => put_line("bad input data, terminating at:"); put_line(Title(1..Title_length)); when finished => put_line("TEST_QR_LINEAR 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_QR_LINEAR;