File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices / ts_qrd.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:24 2000 UTC (23 years, 8 months ago) by maekawa
Branch: PHC, MAIN
CVS Tags: v2, maekawa-ipv6, RELEASE_1_2_3, RELEASE_1_2_2_KNOPPIX_b, RELEASE_1_2_2_KNOPPIX, RELEASE_1_2_2, RELEASE_1_2_1, HEAD Changes since 1.1: +0 -0
lines
Import the second public release of PHCpack.
OKed by Jan Verschelde.
|
with text_io,integer_io; use text_io,integer_io;
with Standard_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
with Standard_Natural_Vectors;
with Standard_Natural_Vectors_io; use Standard_Natural_Vectors_io;
with Standard_Floating_Vectors;
with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
with Standard_Floating_Matrices;
with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
with Standard_Random_Vectors; use Standard_Random_Vectors;
with Standard_Random_Matrices; use Standard_Random_Matrices;
with Standard_Floating_QR_Decomposition; use Standard_Floating_QR_Decomposition;
with Standard_Floating_Least_Squares; use Standard_Floating_Least_Squares;
with Standard_Complex_Vectors;
with Standard_Complex_Vectors_io; use Standard_Complex_Vectors_io;
with Standard_Complex_Matrices;
with Standard_Complex_Matrices_io; use Standard_Complex_Matrices_io;
with Standard_Complex_QR_Decomposition; use Standard_Complex_QR_Decomposition;
with Standard_Complex_Least_Squares; use Standard_Complex_Least_Squares;
procedure ts_qrd is
-- DESCRIPTION :
-- This program tests the implementation of the QR decomposition
-- and least squares approximation.
-- GENERAL TESTS ON QR DECOMPOSITION :
function Extract_Upper_Triangular
( a : Standard_Floating_Matrices.Matrix )
return Standard_Floating_Matrices.Matrix is
-- DESCRIPTION :
-- Returns the upper triangular part of the matrix a.
res : Standard_Floating_Matrices.Matrix(a'range(1),a'range(2));
begin
for i in a'range(1) loop
for j in a'first(2)..(i-1) loop
res(i,j) := 0.0;
end loop;
for j in i..a'last(2) loop
res(i,j) := a(i,j);
end loop;
end loop;
return res;
end Extract_Upper_Triangular;
function Extract_Upper_Triangular
( a : Standard_Complex_Matrices.Matrix )
return Standard_Complex_Matrices.Matrix is
-- DESCRIPTION :
-- Returns the upper triangular part of the matrix a.
res : Standard_Complex_Matrices.Matrix(a'range(1),a'range(2));
begin
for i in a'range(1) loop
for j in a'first(2)..(i-1) loop
res(i,j) := Create(0.0);
end loop;
for j in i..a'last(2) loop
res(i,j) := a(i,j);
end loop;
end loop;
return res;
end Extract_Upper_Triangular;
function Differences ( a,b : in Standard_Floating_Matrices.Matrix )
return double_float is
-- DESCRIPTION :
-- Returns the sum of the differences of all elements |a(i,j)-b(i,j)|.
sum : double_float := 0.0;
begin
for i in a'range(1) loop
for j in a'range(2) loop
sum := sum + abs(a(i,j)-b(i,j));
end loop;
end loop;
return sum;
end Differences;
function Differences ( a,b : in Standard_Complex_Matrices.Matrix )
return double_float is
-- DESCRIPTION :
-- Returns the sum of the differences of all elements |a(i,j)-b(i,j)|.
sum : double_float := 0.0;
begin
for i in a'range(1) loop
for j in a'range(2) loop
sum := sum + AbsVal(a(i,j)-b(i,j));
end loop;
end loop;
return sum;
end Differences;
procedure Orthogonality ( q : in Standard_Floating_Matrices.Matrix ) is
-- DESCRIPTION :
-- Tests whether the columns are orthogonal w.r.t. each other.
sum,ip : double_float;
begin
sum := 0.0;
for j in q'range(2) loop
for k in j+1..q'last(2) loop
ip := 0.0;
for i in q'range(1) loop
ip := ip + q(i,j)*q(i,k);
end loop;
sum := sum + abs(ip);
end loop;
end loop;
put("Orthogonality check : "); put(sum,3,3,3); new_line;
end Orthogonality;
procedure Orthogonality ( q : in Standard_Complex_Matrices.Matrix ) is
-- DESCRIPTION :
-- Tests whether the columns are orthogonal w.r.t. each other.
sum : double_float := 0.0;
ip : Complex_Number;
begin
for j in q'range(2) loop
for k in j+1..q'last(2) loop
ip := Create(0.0);
for i in q'range(1) loop
ip := ip + Conjugate(q(i,j))*q(i,k);
end loop;
sum := sum + AbsVal(ip);
end loop;
end loop;
put("Orthogonality check : "); put(sum,3,3,3); new_line;
end Orthogonality;
procedure Test_QRD ( a,q,r : in Standard_Floating_Matrices.Matrix ) is
wrk : Standard_Floating_Matrices.Matrix(a'range(1),a'range(2));
use Standard_Floating_Matrices;
begin
put_line("The upper triangular part R :"); put(r,3);
wrk := q*r;
put_line("q*r :"); put(wrk,3);
put("Difference in 1-norm between the matrix and q*r: ");
put(Differences(a,wrk),3,3,3); new_line;
Orthogonality(q);
end Test_QRD;
procedure Test_QRD ( a,q,r : in Standard_Complex_Matrices.Matrix ) is
wrk : Standard_Complex_Matrices.Matrix(a'range(1),a'range(2));
use Standard_Complex_Matrices;
begin
put_line("The upper triangular part R :"); put(r,3);
wrk := q*r;
put_line("q*r :"); put(wrk,3);
put("Difference in 1-norm between the matrix and q*r : ");
put(Differences(a,wrk),3,3,3); new_line;
Orthogonality(q);
end Test_QRD;
-- REAL TEST DRIVERS :
procedure Real_LS_Test ( n,m : in natural; piv : in boolean;
a : in Standard_Floating_Matrices.Matrix;
b : in Standard_Floating_Vectors.Vector ) is
wrk : Standard_Floating_Matrices.Matrix(1..n,1..m) := a;
qraux : Standard_Floating_Vectors.Vector(1..m) := (1..m => 0.0);
jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
sol : Standard_Floating_Vectors.Vector(1..m);
rsd,dum : Standard_Floating_Vectors.Vector(1..n);
info : integer;
use Standard_Floating_Matrices;
use Standard_Floating_Vectors;
begin
put_line("The matrix : "); put(a,3);
QRD(wrk,qraux,jpvt,piv);
put_line("The matrix after QR : "); put(wrk,3);
put_line("The vector qraux : "); put(qraux,3); new_line;
if piv
then put("The vector jpvt : "); put(jpvt); new_line;
Permute_Columns(wrk,jpvt);
end if;
QRLS(wrk,n,n,m,qraux,b,dum,dum,sol,rsd,dum,110,info);
if piv
then Permute(sol,jpvt);
end if;
put_line("The solution : "); put(sol,3); new_line;
put_line("The residual : "); put(rsd,3); new_line;
dum := b - a*sol;
put_line("right-hand size - matrix*solution : "); put(dum,3); new_line;
end Real_LS_Test;
procedure Real_QR_Test ( n,m : in natural; piv : in boolean;
a : in Standard_Floating_Matrices.Matrix ) is
wrk : Standard_Floating_Matrices.Matrix(1..n,1..m) := a;
bas : Standard_Floating_Matrices.Matrix(1..n,1..n);
qraux : Standard_Floating_Vectors.Vector(1..m) := (1..m => 0.0);
jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
begin
put_line("The matrix : "); put(a,3);
QRD(wrk,qraux,jpvt,piv);
put_line("The matrix after QR : "); put(wrk,3);
put_line("The vector qraux : "); put(qraux,3); new_line;
if piv
then put("The vector jpvt : "); put(jpvt); new_line;
Permute_Columns(wrk,jpvt);
end if;
for i in wrk'range(1) loop
for j in wrk'range(2) loop
bas(i,j) := wrk(i,j);
end loop;
for j in n+1..m loop
bas(i,j) := 0.0;
end loop;
end loop;
Basis(bas,a);
put_line("The orthogonal part Q of QR :"); put(bas,3);
Test_QRD(a,bas,Extract_Upper_Triangular(wrk));
end Real_QR_Test;
procedure Interactive_Real_QR_Test ( n,m : in natural; piv : in boolean ) is
a : Standard_Floating_Matrices.Matrix(1..n,1..m);
begin
put("Give a "); put(n,1); put("x"); put(m,1);
put_line(" matrix : "); get(a);
Real_QR_Test(n,m,piv,a);
end Interactive_Real_QR_Test;
procedure Interactive_Real_LS_Test ( n,m : in natural; piv : in boolean ) is
a : Standard_Floating_Matrices.Matrix(1..n,1..m);
b : Standard_Floating_Vectors.Vector(1..n);
begin
put("Give a "); put(n,1); put("x"); put(m,1);
put_line(" matrix : "); get(a);
put("Give right-hand size "); put(n,1);
put_line("-vector : "); get(b);
Real_LS_Test(n,m,piv,a,b);
end Interactive_Real_LS_Test;
procedure Random_Real_QR_Test ( n,m : in natural; piv : in boolean ) is
a : Standard_Floating_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
begin
Real_QR_Test(n,m,piv,a);
end Random_Real_QR_Test;
procedure Random_Real_LS_Test ( n,m : in natural; piv : in boolean ) is
a : Standard_Floating_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
b : Standard_Floating_Vectors.Vector(1..n) := Random_Vector(n);
begin
Real_LS_Test(n,m,piv,a,b);
end Random_Real_LS_Test;
-- COMPLEX TEST DRIVERS :
procedure Complex_QR_Test ( n,m : in natural; piv : in boolean;
a : Standard_Complex_Matrices.Matrix ) is
wrk : Standard_Complex_Matrices.Matrix(1..n,1..m) := a;
bas : Standard_Complex_Matrices.Matrix(1..n,1..n);
qraux : Standard_Complex_Vectors.Vector(1..m) := (1..m => Create(0.0));
jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
begin
put_line("The matrix : "); put(a,3);
QRD(wrk,qraux,jpvt,piv);
put_line("The matrix after QR : "); put(wrk,3);
put_line("The vector qraux : "); put(qraux,3); new_line;
-- put("The vector jpvt : "); put(jpvt); new_line;
if not piv
then for i in wrk'range(1) loop
for j in wrk'range(2) loop
bas(i,j) := wrk(i,j);
end loop;
for j in n+1..m loop
bas(i,j) := Create(0.0);
end loop;
end loop;
Basis(bas,a);
put_line("The orthogonal part Q of QR :"); put(bas,3);
Test_QRD(a,bas,Extract_Upper_Triangular(wrk));
end if;
end Complex_QR_Test;
procedure Complex_LS_Test ( n,m : in natural; piv : in boolean;
a : Standard_Complex_Matrices.Matrix;
b : Standard_Complex_Vectors.Vector ) is
wrk : Standard_Complex_Matrices.Matrix(1..n,1..m) := a;
qraux : Standard_Complex_Vectors.Vector(1..m) := (1..m => Create(0.0));
jpvt : Standard_Natural_Vectors.Vector(1..m) := (1..m => 0);
sol : Standard_Complex_Vectors.Vector(1..m);
rsd,dum : Standard_Complex_Vectors.Vector(1..n);
info : integer;
use Standard_Complex_Matrices;
use Standard_Complex_Vectors;
begin
put_line("The matrix : "); put(a,3);
QRD(wrk,qraux,jpvt,piv);
put_line("The matrix after QR : "); put(wrk,3);
put_line("The vector qraux : "); put(qraux,3); new_line;
-- put("The vector jpvt : "); put(jpvt); new_line;
QRLS(wrk,n,n,m,qraux,b,dum,dum,sol,rsd,dum,110,info);
put_line("The solution : "); put(sol,3); new_line;
put_line("The residual : "); put(rsd,3); new_line;
dum := b - a*sol;
put_line("right-hand size - matrix*solution : "); put(dum,3); new_line;
end Complex_LS_Test;
procedure Interactive_Complex_QR_Test
( n,m : in natural; piv : in boolean ) is
a : Standard_Complex_Matrices.Matrix(1..n,1..m);
begin
put("Give a "); put(n,1); put("x"); put(m,1);
put_line(" matrix : "); get(a);
Complex_QR_Test(n,m,piv,a);
end Interactive_Complex_QR_Test;
procedure Interactive_Complex_LS_Test
( n,m : in natural; piv : in boolean ) is
a : Standard_Complex_Matrices.Matrix(1..n,1..m);
b : Standard_Complex_Vectors.Vector(1..n);
begin
put("Give a "); put(n,1); put("x"); put(m,1);
put_line(" matrix : "); get(a);
put("Give right-hand size "); put(n,1);
put_line("-vector : "); get(b);
Complex_LS_Test(n,m,piv,a,b);
end Interactive_Complex_LS_Test;
procedure Random_Complex_QR_Test ( n,m : in natural; piv : in boolean ) is
a : Standard_Complex_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
begin
Complex_QR_Test(n,m,piv,a);
end Random_Complex_QR_Test;
procedure Random_Complex_LS_Test ( n,m : in natural; piv : in boolean ) is
a : Standard_Complex_Matrices.Matrix(1..n,1..m) := Random_Matrix(n,m);
b : Standard_Complex_Vectors.Vector(1..n) := Random_Vector(n);
begin
Complex_LS_Test(n,m,piv,a,b);
end Random_Complex_LS_Test;
-- MAIN PROGRAM :
procedure Main is
n,m : natural;
ans,choice : character;
piv : boolean := false;
begin
new_line;
put_line("Test on the QR decomposition and Least Squares.");
loop
new_line;
put_line("Choose one of the following : ");
put_line(" 0. Exit this program.");
put_line(" 1. QR-decomposition on user-given floating-point matrix.");
put_line(" 2. complex matrix.");
put_line(" 3. on random floating-point matrix.");
put_line(" 4. complex matrix.");
put_line(" 5. Least Squares on user-given floating-point matrix.");
put_line(" 6. complex matrix.");
put_line(" 7. on random floating-point matrix.");
put_line(" 8. complex matrix.");
put("Make your choice (0,1,2,3,4,5,6,7 or 8) : "); get(choice);
exit when (choice = '0');
put("Give the number of rows : "); get(n);
put("Give the number of columns : "); get(m);
put("Do you want pivoting ? (y/n) "); get(ans);
piv := (ans = 'y');
case choice is
when '1' => Interactive_Real_QR_Test(n,m,piv);
when '2' => Interactive_Complex_QR_Test(n,m,piv);
when '3' => Random_Real_QR_Test(n,m,piv);
when '4' => Random_Complex_QR_Test(n,m,piv);
when '5' => Interactive_Real_LS_Test(n,m,piv);
when '6' => Interactive_Complex_LS_Test(n,m,piv);
when '7' => Random_Real_LS_Test(n,m,piv);
when '8' => Random_Complex_LS_Test(n,m,piv);
when others => null;
end case;
end loop;
end Main;
begin
Main;
end ts_qrd;