File: [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices / standard_complex_least_squares.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 Standard_Floating_Numbers; use Standard_Floating_Numbers;
with Standard_Complex_Numbers; use Standard_Complex_Numbers;
package body Standard_Complex_Least_Squares is
-- AUXILIARIES :
function min0 ( a,b : integer ) return integer is
-- DESCRIPTION : returns the minimum of a and b.
begin
if a <= b
then return a;
else return b;
end if;
end min0;
procedure zcopy ( n,start : in natural;
x : in Standard_Complex_Vectors.Vector;
y : out Standard_Complex_Vectors.Vector ) is
-- DESCRIPTION :
-- Copies the n entries of x to y, beginning at start.
begin
for i in start..start+n-1 loop
y(i) := x(i);
end loop;
end zcopy;
function zdotc ( row : natural;
x : Standard_Complex_Matrices.Matrix;
y : Standard_Complex_Vectors.Vector )
return Complex_Number is
-- DESCRIPTION :
-- Dot product of two vectors : x(row..x'last(1),row)*y(row..y'last).
res : Complex_Number := Create(0.0);
begin
for i in row..y'last loop
res := res + Conjugate(x(i,row))*y(i);
end loop;
return res;
end zdotc;
procedure zaxpy ( n,row,col : in natural; f : in Complex_Number;
x : in Standard_Complex_Matrices.Matrix;
y : in out Standard_Complex_Vectors.Vector ) is
-- DESCRIPTION :
-- y(i) := y(i) + f*x(i,col) for i in row..row+n-1.
begin
for i in row..row+n-1 loop
y(i) := y(i) + f*x(i,col);
end loop;
end zaxpy;
procedure QRLS ( x : in out Standard_Complex_Matrices.Matrix;
ldx,n,k : in natural;
qraux,y : in Standard_Complex_Vectors.Vector;
qy,qty,b,rsd,xb : out Standard_Complex_Vectors.Vector;
job : in natural; info : out natural ) is
cb,cqy,cqty,cr,cxb : boolean;
jj,ju,kp1 : integer;
t,temp : Complex_Number;
begin
info := 0; -- set info flag
cqy := (job/10000 /= 0); -- determine what to compute
cqty := (job mod 10000 /= 0);
cb := ((job mod 1000)/100 /= 0);
cr := ((job mod 100)/10 /= 0);
cxb := ((job mod 10) /= 0);
ju := min0(k,n-1);
if ju = 0 -- special action when n = 1
then if cqy then qy(1) := y(1); end if;
if cqty then qty(1) := y(1); end if;
if cxb then xb(1) := y(1); end if;
if cb
then if AbsVal(x(1,1)) = 0.0
then info := 1;
else b(1) := y(1)/x(1,1);
end if;
end if;
if cr then rsd(1) := Create(0.0); end if;
return;
end if;
if cqy -- set up to compute qy or qty
then zcopy(n,y'first,y,qy);
end if;
if cqty
then zcopy(n,y'first,y,qty);
end if;
if cqy -- compute qy
then for j in 1..ju loop
jj := ju - j + 1;
if AbsVal(qraux(jj)) /= 0.0
then temp := x(jj,jj);
x(jj,jj) := qraux(jj);
t := -zdotc(jj,x,qy)/x(jj,jj);
zaxpy(n-jj+1,jj,jj,t,x,qy);
x(jj,jj) := temp;
end if;
end loop;
end if;
if cqty -- compute trans(q)*y
then for j in 1..ju loop
if AbsVal(qraux(j)) /= 0.0
then temp := x(j,j);
x(j,j) := qraux(j);
t := -zdotc(j,x,qty)/x(j,j);
zaxpy(n-j+1,j,j,t,x,qty);
x(j,j) := temp;
end if;
end loop;
end if;
if cb -- set up to compute b,rsd, or xb
then zcopy(k,qty'first,qty,b);
end if;
kp1 := k + 1;
if cxb then zcopy(k,qty'first,qty,xb); end if;
if (cr and (k < n))
then zcopy(n-k,kp1,qty,rsd);
end if;
if (cxb and (kp1 <= n))
then for i in kp1..n loop
xb(i) := Create(0.0);
end loop;
end if;
if cr
then for i in 1..k loop
rsd(i) := Create(0.0);
end loop;
end if;
if cb -- compute b
then for j in 1..k loop
jj := k - j + 1;
if AbsVal(x(jj,jj)) = 0.0
then info := jj; exit;
end if;
b(jj) := b(jj)/x(jj,jj);
if jj /= 1
then t := -b(jj);
zaxpy(jj-1,1,jj,t,x,b);
end if;
end loop;
end if;
if cr or cxb -- compute rsd or xb as requested
then for j in 1..ju loop
jj := ju - j + 1;
if AbsVal(qraux(jj)) /= 0.0
then temp := x(jj,jj);
x(jj,jj) := qraux(jj);
if cr
then t := -zdotc(jj,x,rsd)/x(jj,jj);
zaxpy(n-jj+1,jj,jj,t,x,rsd);
end if;
if cxb
then t := -zdotc(jj,x,xb)/x(jj,jj);
zaxpy(n-jj+1,jj,jj,t,x,xb);
end if;
x(jj,jj) := temp;
end if;
end loop;
end if;
end QRLS;
end Standard_Complex_Least_Squares;