File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices / multprec_complex_linear_solvers.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 Multprec_Complex_Numbers; use Multprec_Complex_Numbers;
package body Multprec_Complex_Linear_Solvers is
-- AUXLILIARIES :
function dconjg ( x : Complex_Number ) return Complex_Number is
-- DESCRIPTION :
-- Returns the complex conjugate of x.
res : Complex_Number;
re : Floating_Number := REAL_PART(x);
im : Floating_Number := IMAG_PART(x);
begin
Min(im);
res := Create(re,im);
Clear(re); Clear(im);
return res;
end dconjg;
function csign ( x,y : Complex_Number ) return Complex_Number is
-- DESCRIPTION :
-- Returns |x|*y/|y|.
res : Complex_Number;
fltacc : Floating_Number;
cmpacc : Complex_Number;
begin
fltacc := AbsVal(x);
res := Create(fltacc);
Clear(fltacc);
Mul(res,y);
fltacc := AbsVal(y);
cmpacc := Create(fltacc);
Div(res,cmpacc);
Clear(fltacc);
Clear(cmpacc);
return res;
end csign;
function Inverse_Abs_Sum ( z : Vector ) return Floating_Number is
-- DESCRIPTION :
-- Returns the reciprocal of the sum of the absolute values in z.
res,sum,acc : Floating_Number;
begin
sum := Create(0);
for i in z'range loop
acc := AbsVal(z(i));
Add(sum,acc);
Clear(acc);
end loop;
res := Create(1);
Div(res,sum);
Clear(sum);
return res;
end Inverse_Abs_Sum;
-- 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 : Floating_Number := AbsVal(a(i,res));
tmp : Floating_Number;
begin
for j in a'first(2)+1..a'last(2) loop
tmp := AbsVal(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,fltacc : Floating_Number;
temp,cmpacc : Complex_Number;
begin
info := 0;
nm1 := n - 1;
if nm1 >= 1
then for k in 1..nm1 loop
kp1 := k + 1;
l := k;
smax := AbsVal(a(k,k)); -- find the pivot index l
for i in kp1..n loop
fltacc := AbsVal(a(i,k));
if fltacc > smax
then l := i;
Copy(fltacc,smax);
end if;
Clear(fltacc);
end loop;
ipvt(k) := l;
if Equal(smax,0.0) -- this column is already triangularized
then info := k;
else if l /= k -- interchange if necessary
then Copy(a(l,k),temp);
Copy(a(k,k),a(l,k));
Copy(temp,a(k,k)); Clear(temp);
end if;
cmpacc := Create(1); -- compute multipliers
Div(cmpacc,a(k,k));
Min(cmpacc);
for i in kp1..n loop
Mul(a(i,k),cmpacc);
end loop;
Clear(cmpacc);
for j in kp1..n loop -- row elimination
Copy(a(l,j),temp);
if l /= k
then Copy(a(k,j),a(l,j));
Copy(temp,a(k,j));
end if;
for i in kp1..n loop
cmpacc := temp*a(i,k);
Add(a(i,j),cmpacc);
Clear(cmpacc);
end loop;
end loop;
Clear(temp);
end if;
Clear(smax);
end loop;
end if;
ipvt(n) := n;
fltacc := AbsVal(a(n,n));
if Equal(fltacc,0.0)
then info := n;
end if;
Clear(fltacc);
end lufac;
procedure lufco ( a : in out Matrix; n : in integer;
ipvt : out Standard_Natural_Vectors.Vector;
rcond : out Floating_Number ) 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 : Multprec_Complex_Vectors.Vector(1..n);
info,kb,kp1,l : integer;
s,sm,sum,anorm,ynorm,fltacc1,fltacc2 : Floating_Number;
ek,t,wk,wkm,cmpacc1,cmpacc2 : Complex_Number;
ipvtt : Standard_Natural_Vectors.Vector(1..n);
begin
anorm := Create(0); -- compute 1-norm of a
for j in 1..n loop
sum := Create(0);
for i in 1..n loop
fltacc1 := AbsVal(a(i,j));
Add(sum,fltacc1);
Clear(fltacc1);
end loop;
if sum > anorm
then Copy(sum,anorm);
end if;
Clear(sum);
end loop;
lufac(a,n,ipvtt,info); -- factor
ipvt := ipvtt;
ek := Create(1); -- solve ctrans(u)*w = e
for j in 1..n loop
z(j) := Create(0);
end loop;
for k in 1..n loop
fltacc1 := AbsVal(z(k));
if not Equal(fltacc1,0.0)
then cmpacc1 := -z(k);
Clear(ek);
ek := csign(ek,cmpacc1);
Clear(cmpacc1);
end if;
Clear(fltacc1);
cmpacc1 := ek - z(k);
fltacc1 := AbsVal(cmpacc1);
Clear(cmpacc1);
fltacc2 := AbsVal(a(k,k));
if fltacc1 > fltacc2
then s := fltacc2/fltacc1;
cmpacc1 := Create(s);
Mul(z,cmpacc1);
Mul(ek,cmpacc1);
Clear(cmpacc1);
Clear(s);
end if;
Clear(fltacc1); Clear(fltacc2);
Clear(wk); wk := ek - z(k);
Clear(wkm); wkm := ek + z(k);
Min(wkm);
s := AbsVal(wk);
sm := AbsVal(wkm);
fltacc1 := AbsVal(a(k,k));
if Equal(fltacc1,0.0)
then Clear(wk); wk := Create(1);
Clear(wkm); wkm := Create(1);
else cmpacc1 := dconjg(a(k,k));
Div(wk,cmpacc1);
Div(wkm,cmpacc1);
Clear(cmpacc1);
end if;
Clear(fltacc1);
kp1 := k + 1;
if kp1 <= n
then for j in kp1..n loop
cmpacc2 := dconjg(a(k,j));
cmpacc1 := wkm*cmpacc2;
Add(cmpacc1,z(j));
fltacc1 := AbsVal(cmpacc1);
Add(sm,fltacc1);
Clear(fltacc1);
Clear(cmpacc1);
cmpacc1 := wk*cmpacc2;
Add(z(j),cmpacc1);
Clear(cmpacc1);
Clear(cmpacc2);
fltacc1 := AbsVal(z(j));
Add(s,fltacc1);
Clear(fltacc1);
end loop;
if s < sm
then t := wkm - wk;
Copy(wkm,wk);
for j in kp1..n loop
cmpacc1 := t*dconjg(a(k,j));
Add(z(j),cmpacc1);
Clear(cmpacc1);
end loop;
Clear(t);
end if;
end if;
Copy(wk,z(k));
end loop;
s := Inverse_Abs_Sum(z);
cmpacc1 := Create(s);
Mul(z,cmpacc1);
Clear(cmpacc1);
Clear(s);
for k in 1..n loop -- solve ctrans(l)*y = w
kb := n+1-k;
if kb < n
then t := Create(0);
for i in (kb+1)..n loop
cmpacc1 := dconjg(a(i,kb))*z(i);
Add(t,cmpacc1);
Clear(cmpacc1);
end loop;
Add(z(kb),t);
end if;
fltacc1 := AbsVal(z(kb));
if fltacc1 > 1.0
then s := Create(1);
Div(s,fltacc1);
cmpacc1 := Create(s);
Mul(z,cmpacc1);
Clear(cmpacc1);
Clear(s);
end if;
Clear(fltacc1);
l := ipvtt(kb);
if l /= kb
then Copy(z(l),t);
Copy(z(kb),z(l));
Copy(t,z(kb));
Clear(t);
end if;
end loop;
s := Inverse_Abs_Sum(z);
cmpacc1 := Create(s);
Clear(s);
Mul(z,cmpacc1);
Clear(cmpacc1);
ynorm := Create(1);
for k in 1..n loop -- solve l*v = y
l := ipvtt(k);
if l /= k
then Copy(z(l),t);
Copy(z(k),z(l));
Copy(t,z(k));
else Copy(z(l),t);
end if;
if k < n
then for i in (k+1)..n loop
cmpacc1 := t*a(i,k);
Add(z(i),cmpacc1);
Clear(cmpacc1);
end loop;
end if;
Clear(t);
fltacc1 := AbsVal(z(k));
if fltacc1 > 1.0
then s := Create(1);
Div(s,fltacc1);
cmpacc1 := Create(s);
Mul(z,cmpacc1);
Clear(cmpacc1);
Mul(ynorm,s);
Clear(s);
end if;
Clear(fltacc1);
end loop;
s := Inverse_Abs_Sum(z);
cmpacc1 := Create(s);
Mul(z,cmpacc1);
Clear(cmpacc1);
Mul(ynorm,s);
Clear(s);
for k in 1..n loop -- solve u*z = v
kb := n+1-k;
fltacc1 := AbsVal(z(kb));
fltacc2 := AbsVal(a(kb,kb));
if fltacc1 > fltacc2
then s := fltacc2/fltacc1;
cmpacc1 := Create(s);
Mul(z,cmpacc1);
Mul(ynorm,s);
Clear(s);
end if;
if Equal(fltacc2,0.0)
then Clear(z(kb));
z(kb) := Create(1);
else Div(z(kb),a(kb,kb));
end if;
Clear(fltacc1);
Clear(fltacc2);
t := -z(kb);
for i in 1..(kb-1) loop
cmpacc1 := t*a(i,kb);
Add(z(i),cmpacc1);
Clear(cmpacc1);
end loop;
Clear(t);
end loop;
s := Inverse_Abs_Sum(z); -- make znorm = 1.0
cmpacc1 := Create(s);
Mul(z,cmpacc1);
Clear(cmpacc1);
Mul(ynorm,s);
if Equal(anorm,0.0)
then rcond := Create(0);
else rcond := ynorm/anorm;
end if;
Clear(ynorm); Clear(s); Clear(sm);
Clear(ek); Clear(t); Clear(wk); Clear(wkm);
Clear(z);
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,acc : Complex_Number;
begin
nm1 := n-1;
if nm1 >= 1 -- solve l*y = b
then for k in 1..nm1 loop
l := ipvt(k);
Copy(b(l),temp);
if l /= k
then Copy(b(k),b(l));
Copy(temp,b(k));
end if;
for i in (k+1)..n loop
acc := temp*a(i,k);
Add(b(i),acc);
Clear(acc);
end loop;
Clear(temp);
end loop;
end if;
for k in 1..n loop -- solve u*x = y
kb := n+1-k;
Div(b(kb),a(kb,kb));
temp := -b(kb);
for j in 1..(kb-1) loop
acc := temp*a(j,kb);
Add(b(j),acc);
Clear(acc);
end loop;
Clear(temp);
end loop;
end lusolve;
procedure Triangulate ( a : in out Matrix; n,m : in integer ) is
max,cbs : Floating_Number;
temp : Complex_Number;
pivot,k,kcolumn : integer;
begin
k := 1;
kcolumn := 1;
while (k <= n) and (kcolumn <= m) loop
max := Create(0); -- find pivot
pivot := 0;
for l in k..n loop
cbs := AbsVal(a(l,kcolumn));
if 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);
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);
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 : Floating_Number;
temp : Complex_Number;
pivot,k,kcolumn : integer;
begin
k := 1;
kcolumn := 1;
while (k <= n) and (kcolumn <= m) loop
max := Create(0); -- find pivot
for l in k..n loop
if AbsVal(a(l,kcolumn)) > max
then max := AbsVal(a(l,kcolumn));
pivot := l;
end if;
end loop;
if Equal(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);
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);
end loop;
for j in (k+1)..n loop
a(j,kcolumn) := Create(0);
end loop;
k := k + 1;
kcolumn := kcolumn + 1;
end if;
end loop;
end Diagonalize;
end Multprec_Complex_Linear_Solvers;