File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports / givens_rotations.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:27 2000 UTC (23 years, 10 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 Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
package body Givens_Rotations is
procedure Givens_Factors ( v: in Vector; i,j : in integer;
cosi,sine : out double_float ) is
norm : double_float;
begin
norm := SQRT(v(i)*v(i) + v(j)*v(j));
cosi := v(i)/norm; sine := v(j)/norm;
end Givens_Factors;
procedure Givens_Factors ( mat : in Matrix; i,j : in integer;
cosi,sine : out double_float ) is
norm : double_float;
begin
norm := SQRT(mat(i,i)*mat(i,i) + mat(j,i)*mat(j,i));
cosi := mat(i,i)/norm; sine := mat(j,i)/norm;
end Givens_Factors;
procedure Givens_Rotation ( v : in out Vector; i,j : in integer;
cosi,sine : in double_float ) is
temp : double_float;
begin
temp := v(i);
v(i) := cosi*v(i) + sine*v(j);
v(j) := -sine*temp + cosi*v(j);
end Givens_Rotation;
procedure Givens_Rotation ( mat : in out Matrix; i,j : in integer;
cosi,sine : in double_float ) is
temp : double_float;
begin
-- for k in mat'range(2) loop -- certain angle preservation
for k in i..mat'last(2) loop -- only angle preservation if upper triangular
temp := mat(i,k);
mat(i,k) := cosi*mat(i,k) + sine*mat(j,k);
mat(j,k) := -sine*temp + cosi*mat(j,k);
end loop;
end Givens_Rotation;
procedure Givens_Rotation ( mat : in out Matrix; lastcol,i,j : in integer;
cosi,sine : in double_float ) is
temp : double_float;
begin
-- for k in mat'first(2)..lastcol loop -- certain angle preservation
for k in i..lastcol loop -- only angle preservation if upper triangular
temp := mat(i,k);
mat(i,k) := cosi*mat(i,k) + sine*mat(j,k);
mat(j,k) := -sine*temp + cosi*mat(j,k);
end loop;
end Givens_Rotation;
procedure Givens_Rotation ( v : in out Vector; i,j : in integer ) is
cosi,sine : double_float;
begin
Givens_Factors(v,i,j,cosi,sine);
Givens_Rotation(v,i,j,cosi,sine);
end Givens_Rotation;
procedure Givens_Rotation ( mat : in out Matrix; i,j : in integer ) is
cosi,sine : double_float;
begin
Givens_Factors(mat,i,j,cosi,sine);
Givens_Rotation(mat,i,j,cosi,sine);
end Givens_Rotation;
procedure Givens_Rotation ( mat : in out Matrix;
lastcol,i,j : in integer ) is
cosi,sine : double_float;
begin
Givens_Factors(mat,i,j,cosi,sine);
Givens_Rotation(mat,lastcol,i,j,cosi,sine);
end Givens_Rotation;
procedure Upper_Triangulate
( row : in integer; mat : in out Matrix; tol : in double_float;
ipvt : in out Standard_Integer_Vectors.Vector;
pivot : out integer ) is
piv,tpi : integer := 0;
tmp : double_float;
begin
for j in row..mat'last(2) loop -- search pivot
if abs(mat(row,j)) > tol
then piv := j;
end if;
exit when (piv /= 0);
end loop;
if piv /= 0 -- zero row
then
if piv /= row -- interchange columns
then for k in mat'range(1) loop
tmp := mat(k,row); mat(k,row) := mat(k,piv); mat(k,piv) := tmp;
end loop;
tpi := ipvt(row); ipvt(row) := ipvt(piv); ipvt(piv) := tpi;
end if;
for k in row+1..mat'last(1) loop -- perform Givens rotations
if abs(mat(k,row)) > tol
then Givens_Rotation(mat,row,k);
end if;
end loop;
end if;
pivot := piv;
end Upper_Triangulate;
procedure Upper_Triangulate ( mat : in out Matrix; tol : in double_float ) is
pivot : integer;
temp : double_float;
begin
for i in mat'range(1) loop -- mat(mat'first(1)..i-1,mat'first(2)..i-1)
pivot := 0; -- is already upper triangular
for j in i..mat'last(2) loop -- search pivot
if abs(mat(i,j)) > tol
then pivot := j;
end if;
exit when (pivot /= 0);
end loop;
exit when (pivot = 0); -- zero row
if pivot /= i -- interchange columns
then for k in mat'range(1) loop
temp := mat(k,i); mat(k,i) := mat(k,pivot); mat(k,pivot) := temp;
end loop;
end if;
for k in i+1..mat'last(1) loop -- perform Givens rotations
if abs(mat(k,i)) > tol
then Givens_Rotation(mat,i,k);
end if;
end loop; -- mat(mat'first(1)..i,mat'first(2)..i)
exit when i > mat'last(2); -- is upper triangular
end loop;
end Upper_Triangulate;
procedure Upper_Triangulate ( mat : in out Matrix; tol : in double_float;
ipvt : out Standard_Integer_Vectors.Vector ) is
pivot,tpi : integer;
pivots : Standard_Integer_Vectors.Vector(mat'range(2));
temp : double_float;
begin
for i in pivots'range loop -- initialize vector of pivots
pivots(i) := i;
end loop;
for i in mat'range(1) loop -- mat(mat'first(1)..i-1,mat'first(2)..i-1)
pivot := 0; -- is already upper triangular
for j in i..mat'last(2) loop -- search pivot
if abs(mat(i,j)) > tol
then pivot := j;
end if;
exit when (pivot /= 0);
end loop;
exit when (pivot = 0); -- zero row
if pivot /= i -- interchange columns
then for k in mat'range(1) loop
temp := mat(k,i); mat(k,i) := mat(k,pivot); mat(k,pivot) := temp;
end loop;
tpi := pivots(i); pivots(i) := pivots(pivot); pivots(pivot) := tpi;
end if;
for k in i+1..mat'last(1) loop -- perform Givens rotations
if abs(mat(k,i)) > tol
then Givens_Rotation(mat,i,k);
end if;
end loop; -- mat(mat'first(1)..i,mat'first(2)..i)
exit when i > mat'last(2); -- is upper triangular
end loop;
ipvt := pivots;
end Upper_Triangulate;
procedure Upper_Triangulate ( mat : in out Matrix; rhs : in out Vector;
tol : in double_float ) is
pivot : integer;
temp,cosi,sine : double_float;
begin
for i in mat'range(1) loop -- mat(mat'first(1)..i-1,mat'first(2)..i-1)
pivot := 0; -- is already upper triangular
for j in i..mat'last(2) loop -- search pivot
if abs(mat(i,j)) > tol
then pivot := j;
end if;
exit when (pivot /= 0);
end loop;
exit when (pivot = 0); -- zero row
if pivot /= i -- interchange columns
then for k in mat'range(1) loop
temp := mat(k,i); mat(k,i) := mat(k,pivot); mat(k,pivot) := temp;
end loop;
end if;
for k in i+1..mat'last(1) loop -- perform Givens rotations
if abs(mat(k,i)) > tol
then Givens_Factors(mat,i,k,cosi,sine);
Givens_Rotation(mat,i,k,cosi,sine);
Givens_Rotation(rhs,i,k,cosi,sine);
end if;
end loop; -- mat(mat'first(1)..i,mat'first(2)..i)
exit when i > mat'last(2); -- is upper triangular
end loop;
end Upper_Triangulate;
procedure Upper_Triangulate ( mat : in out Matrix; rhs : in out Vector;
tol : in double_float;
ipvt : out Standard_Integer_Vectors.Vector ) is
pivot,tpi : integer;
pivots : Standard_Integer_Vectors.Vector(mat'range(2));
temp,cosi,sine : double_float;
begin
for i in pivots'range loop -- initialize vector of pivots
pivots(i) := i;
end loop;
for i in mat'range(1) loop -- mat(mat'first(1)..i-1,mat'first(2)..i-1)
pivot := 0; -- is already upper triangular
for j in i..mat'last(2) loop -- search pivot
if abs(mat(i,j)) > tol
then pivot := j;
end if;
exit when (pivot /= 0);
end loop;
exit when (pivot = 0); -- zero row
if pivot /= i -- interchange columns
then for k in mat'range(1) loop
temp := mat(k,i); mat(k,i) := mat(k,pivot); mat(k,pivot) := temp;
end loop;
tpi := pivots(i); pivots(i) := pivots(pivot); pivots(pivot) := tpi;
end if;
for k in i+1..mat'last(1) loop -- perform Givens rotations
if abs(mat(k,i)) > tol
then Givens_Factors(mat,i,k,cosi,sine);
Givens_Rotation(mat,i,k,cosi,sine);
Givens_Rotation(rhs,i,k,cosi,sine);
end if;
end loop; -- mat(mat'first(1)..i,mat'first(2)..i)
exit when i > mat'last(2); -- is upper triangular
end loop;
ipvt := pivots;
end Upper_Triangulate;
procedure Solve ( mat : in Matrix; rhs : in Vector; tol : in double_float;
x : out Vector ) is
rank : natural := 0;
sol : vector(mat'range(1)) := (mat'range(1) => 0.0);
begin
for i in mat'range(1) loop -- determine the rank of the matrix
if abs(mat(i,i)) > tol
then rank := rank + 1;
end if;
exit when ((i >= mat'last(2)) or (abs(mat(i,i)) < tol));
end loop;
for i in reverse mat'first(1)..rank loop -- back substitution
for j in i+1..rank loop
sol(i) := sol(i) + mat(i,j)*sol(j);
end loop;
sol(i) := rhs(i) - sol(i);
if abs(mat(i,i)) > tol
then sol(i) := sol(i)/mat(i,i);
elsif abs(sol(i)) < tol
then sol(i) := 1.0;
else return;
end if;
end loop; -- invariant : sol(i..rank) computed
x := sol;
end Solve;
procedure Solve ( mat : in Matrix; rhs : in Vector; tol : in double_float;
rank : in natural; x : out Vector ) is
sol : vector(mat'range(1)) := (mat'range(1) => 0.0);
begin
for i in reverse mat'first(1)..rank loop -- back substitution
for j in i+1..rank loop
sol(i) := sol(i) + mat(i,j)*sol(j);
end loop;
sol(i) := rhs(i) - sol(i);
if abs(mat(i,i)) > tol
then sol(i) := sol(i)/mat(i,i);
elsif abs(sol(i)) < tol
then sol(i) := 1.0;
else return;
end if;
end loop; -- invariant : sol(i..rank) computed
x := sol;
end Solve;
end Givens_Rotations;