File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices / standard_complex_linear_solvers.adb (download)
Revision 1.1.1.1 (vendor branch), Sun Oct 29 17:45:24 2000 UTC (23 years, 9 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_Complex_Numbers; use Standard_Complex_Numbers;
package body Standard_Complex_Linear_Solvers is
-- AUXLILIARIES :
function cabs ( c : Complex_Number ) return double_float is
begin
return (ABS(REAL_PART(c)) + ABS(IMAG_PART(c)));
end cabs;
function dconjg ( x : Complex_Number ) return Complex_Number is
begin
return Create(REAL_PART(x),-IMAG_PART(x));
end dconjg;
function csign ( x,y : Complex_Number ) return Complex_Number is
begin
return (Create(cabs(x)) * y / Create(cabs(y)));
end csign;
-- TARGET ROUTINES :
procedure Scale ( a : in out Matrix; b : in out Vector ) is
fac : Complex_Number;
function Maximum ( a : in Matrix; i : in integer ) return Complex_Number is
res : integer := a'first(2);
max : double_float := cabs(a(i,res));
tmp : double_float;
begin
for j in a'first(2)+1..a'last(2) loop
tmp := cabs(a(i,j));
if tmp > max
then max := tmp; res := j;
end if;
end loop;
return a(i,res);
end Maximum;
procedure Divide ( a : in out Matrix; b : in out Vector;
i : in integer; fac : in Complex_Number ) is
begin
for j in a'range(2) loop
a(i,j) := a(i,j)/fac;
end loop;
b(i) := b(i)/fac;
end Divide;
begin
for i in a'range(1) loop
fac := Maximum(a,i);
Divide(a,b,i,fac);
end loop;
end Scale;
-- TARGET ROUTINES :
procedure lufac ( a : in out Matrix; n : in integer;
ipvt : out Standard_Natural_Vectors.Vector;
info : out natural ) is
kp1,l,nm1 : integer;
smax : double_float;
temp : Complex_Number;
begin
info := 0;
nm1 := n - 1;
if nm1 >= 1
then for k in 1..nm1 loop
kp1 := k + 1;
-- find the pivot index l
l := k; smax := cabs(a(k,k)); --modulus(a(k,k));
for i in kp1..n loop
if cabs(a(i,k)) > smax --modulus(a(i,k)) > smax
then l := i;
smax := cabs(a(i,k)); --modulus(a(i,k));
end if;
end loop;
ipvt(k) := l;
if smax = 0.0
then -- this column is already triangularized
info := k;
else
-- interchange if necessary
if l /= k
then temp := a(l,k);
a(l,k) := a(k,k);
a(k,k) := temp;
end if;
-- compute multipliers
temp := -Create(1.0)/a(k,k);
for i in kp1..n loop
a(i,k) := temp * a(i,k);
end loop;
-- row elimination with column indexing
for j in kp1..n loop
temp := a(l,j);
if l /= k
then a(l,j) := a(k,j);
a(k,j) := temp;
end if;
for i in kp1..n loop
a(i,j) := a(i,j) + temp * a(i,k);
end loop;
end loop;
end if;
end loop;
end if;
ipvt(n) := n;
if AbsVal(a(n,n)) = 0.0
then info := n;
end if;
end lufac;
procedure lufco ( a : in out Matrix; n : in integer;
ipvt : out Standard_Natural_Vectors.Vector;
rcond : out double_float ) is
-- NOTE :
-- rcond = 1/(norm(a)*(estimate of norm(inverse(a))))
-- estimate = norm(z)/norm(y) where a*z = y and ctrans(a)*y = e.
-- ctrans(a) is the conjugate transpose of a.
-- The components of e are chosen to cause maximum local
-- growth in teh elements of w where ctrans(u)*w = e.
-- The vectors are frequently rescaled to avoid overflow.
z : Standard_Complex_Vectors.Vector(1..n);
info,kb,kp1,l : integer;
s,sm,sum,anorm,ynorm : double_float;
ek,t,wk,wkm : Complex_Number;
ipvtt : Standard_Natural_Vectors.Vector(1..n);
begin
anorm := 0.0; -- compute 1-norm of a
for j in 1..n loop
sum := 0.0;
for i in 1..n loop
sum := sum + cabs(a(i,j));
end loop;
if sum > anorm
then anorm := sum;
end if;
end loop;
lufac(a,n,ipvtt,info); -- factor
for i in 1..n loop
ipvt(i) := ipvtt(i);
end loop;
ek := Create(1.0); -- solve ctrans(u)*w = e
for j in 1..n loop
z(j) := Create(0.0);
end loop;
for k in 1..n loop
if cabs(z(k)) /= 0.0
then ek := csign(ek,-z(k));
end if;
if cabs(ek-z(k)) > cabs(a(k,k))
then s := cabs(a(k,k))/cabs(ek-z(k));
z := Create(s) * z;
ek := Create(s) * ek;
end if;
wk := ek - z(k);
wkm := -ek - z(k);
s := cabs(wk);
sm := cabs(wkm);
if cabs(a(k,k)) = 0.0
then wk := Create(1.0);
wkm := Create(1.0);
else wk := wk / dconjg(a(k,k));
wkm := wkm / dconjg(a(k,k));
end if;
kp1 := k + 1;
if kp1 <= n
then for j in kp1..n loop
sm := sm + cabs(z(j)+wkm*dconjg(a(k,j)));
z(j) := z(j) + wk*dconjg(a(k,j));
s := s + cabs(z(j));
end loop;
if s < sm
then t := wkm - wk;
wk := wkm;
for j in kp1..n loop
z(j) := z(j) + t*dconjg(a(k,j));
end loop;
end if;
end if;
z(k) := wk;
end loop;
sum := 0.0;
for i in 1..n loop
sum := sum + cabs(z(i));
end loop;
s := 1.0 / sum;
z := Create(s) * z;
for k in 1..n loop -- solve ctrans(l)*y = w
kb := n+1-k;
if kb < n
then t := Create(0.0);
for i in (kb+1)..n loop
t := t + dconjg(a(i,kb))*z(i);
end loop;
z(kb) := z(kb) + t;
end if;
if cabs(z(kb)) > 1.0
then s := 1.0 / cabs(z(kb));
z := Create(s) * z;
end if;
l := ipvtt(kb);
t := z(l);
z(l) := z(kb);
z(kb) := t;
end loop;
sum := 0.0;
for i in 1..n loop
sum := sum + cabs(z(i));
end loop;
s := 1.0 / sum;
z := Create(s) * z;
ynorm := 1.0;
for k in 1..n loop -- solve l*v = y
l := ipvtt(k);
t := z(l);
z(l) := z(k);
z(k) := t;
if k < n
then for i in (k+1)..n loop
z(i) := z(i) + t * a(i,k);
end loop;
end if;
if cabs(z(k)) > 1.0
then s := 1.0 / cabs(z(k));
z := Create(s) * z;
ynorm := s * ynorm;
end if;
end loop;
sum := 0.0;
for i in 1..n loop
sum := sum + cabs(z(i));
end loop;
s := 1.0 / sum;
z := Create(s) * z;
ynorm := s * ynorm;
for k in 1..n loop -- solve u*z = v
kb := n+1-k;
if cabs(z(kb)) > cabs(a(kb,kb))
then s := cabs(a(kb,kb)) / cabs(z(kb));
z := Create(s) * z;
ynorm := s * ynorm;
end if;
if cabs(a(kb,kb)) = 0.0
then z(kb) := Create(1.0);
else z(kb) := z(kb) / a(kb,kb);
end if;
t := -z(kb);
for i in 1..(kb-1) loop
z(i) := z(i) + t * a(i,kb);
end loop;
end loop;
sum := 0.0; -- make znorm = 1.0
for i in 1..n loop
sum := sum + cabs(z(i));
end loop;
s := 1.0 / sum;
z := Create(s) * z;
ynorm := s * ynorm;
if anorm = 0.0
then rcond := 0.0;
else rcond := ynorm/anorm;
end if;
end lufco;
procedure lusolve ( a : in Matrix; n : in integer;
ipvt : in Standard_Natural_Vectors.Vector;
b : in out Vector ) is
l,nm1,kb : integer;
temp : Complex_Number;
begin
nm1 := n-1;
if nm1 >= 1 -- solve l*y = b
then for k in 1..nm1 loop
l := ipvt(k);
temp := b(l);
if l /= k
then b(l) := b(k);
b(k) := temp;
end if;
for i in (k+1)..n loop
b(i) := b(i) + temp * a(i,k);
end loop;
end loop;
end if;
for k in 1..n loop -- solve u*x = y
kb := n+1-k;
b(kb) := b(kb) / a(kb,kb);
temp := -b(kb);
for j in 1..(kb-1) loop
b(j) := b(j) + temp * a(j,kb);
end loop;
end loop;
end lusolve;
procedure Triangulate ( a : in out Matrix; n,m : in integer ) is
max,cbs : double_float;
temp : Complex_Number;
pivot,k,kcolumn : integer;
tol : constant double_float := 10.0**(-10);
begin
k := 1;
kcolumn := 1;
while (k <= n) and (kcolumn <= m) loop
max := 0.0; -- find pivot
pivot := 0;
for l in k..n loop
cbs := cabs(a(l,kcolumn));
if (cbs > tol) and then (cbs > max)
then max := cbs;
pivot := l;
end if;
end loop;
if pivot = 0
then kcolumn := kcolumn + 1;
else if pivot /= k -- interchange if necessary
then for i in 1..m loop
temp := a(pivot,i);
a(pivot,i) := a(k,i);
a(k,i) := temp;
end loop;
end if;
for j in (kcolumn+1)..m loop -- triangulate a
a(k,j) := a(k,j) / a(k,kcolumn);
end loop;
a(k,kcolumn) := Create(1.0);
for i in (k+1)..n loop
for j in (kcolumn+1)..m loop
a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
end loop;
a(i,kcolumn) := Create(0.0);
end loop;
k := k + 1;
kcolumn := kcolumn + 1;
end if;
end loop;
end Triangulate;
procedure Diagonalize ( a : in out Matrix; n,m : in integer ) is
max : double_float;
temp : Complex_Number;
pivot,k,kcolumn : integer;
begin
k := 1;
kcolumn := 1;
while (k <= n) and (kcolumn <= m) loop
max := 0.0; -- find pivot
for l in k..n loop
if cabs(a(l,kcolumn)) > max
then max := cabs(a(l,kcolumn));
pivot := l;
end if;
end loop;
if max = 0.0
then kcolumn := kcolumn + 1;
else if pivot /= k -- interchange if necessary
then for i in 1..m loop
temp := a(pivot,i);
a(pivot,i) := a(k,i);
a(k,i) := temp;
end loop;
end if;
for j in (kcolumn+1)..m loop -- diagonalize a
a(k,j) := a(k,j) / a(k,kcolumn);
end loop;
a(k,kcolumn) := Create(1.0);
for i in 1..(k-1) loop
for j in (kcolumn+1)..m loop
a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
end loop;
end loop;
for i in (k+1)..n loop
for j in (kcolumn+1)..m loop
a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
end loop;
end loop;
for j in 1..(k-1) loop
a(j,kcolumn) := Create(0.0);
end loop;
for j in (k+1)..n loop
a(j,kcolumn) := Create(0.0);
end loop;
k := k + 1;
kcolumn := kcolumn + 1;
end if;
end loop;
end Diagonalize;
end Standard_Complex_Linear_Solvers;