![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to standard_complex_least_squares.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices|
Return to standard_complex_least_squares.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices |
1.1 maekawa 1: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
2: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
3:
4: package body Standard_Complex_Least_Squares is
5:
6: -- AUXILIARIES :
7:
8: function min0 ( a,b : integer ) return integer is
9:
10: -- DESCRIPTION : returns the minimum of a and b.
11:
12: begin
13: if a <= b
14: then return a;
15: else return b;
16: end if;
17: end min0;
18:
19: procedure zcopy ( n,start : in natural;
20: x : in Standard_Complex_Vectors.Vector;
21: y : out Standard_Complex_Vectors.Vector ) is
22:
23: -- DESCRIPTION :
24: -- Copies the n entries of x to y, beginning at start.
25:
26: begin
27: for i in start..start+n-1 loop
28: y(i) := x(i);
29: end loop;
30: end zcopy;
31:
32: function zdotc ( row : natural;
33: x : Standard_Complex_Matrices.Matrix;
34: y : Standard_Complex_Vectors.Vector )
35: return Complex_Number is
36:
37: -- DESCRIPTION :
38: -- Dot product of two vectors : x(row..x'last(1),row)*y(row..y'last).
39:
40: res : Complex_Number := Create(0.0);
41:
42: begin
43: for i in row..y'last loop
44: res := res + Conjugate(x(i,row))*y(i);
45: end loop;
46: return res;
47: end zdotc;
48:
49: procedure zaxpy ( n,row,col : in natural; f : in Complex_Number;
50: x : in Standard_Complex_Matrices.Matrix;
51: y : in out Standard_Complex_Vectors.Vector ) is
52:
53: -- DESCRIPTION :
54: -- y(i) := y(i) + f*x(i,col) for i in row..row+n-1.
55:
56: begin
57: for i in row..row+n-1 loop
58: y(i) := y(i) + f*x(i,col);
59: end loop;
60: end zaxpy;
61:
62: procedure QRLS ( x : in out Standard_Complex_Matrices.Matrix;
63: ldx,n,k : in natural;
64: qraux,y : in Standard_Complex_Vectors.Vector;
65: qy,qty,b,rsd,xb : out Standard_Complex_Vectors.Vector;
66: job : in natural; info : out natural ) is
67:
68: cb,cqy,cqty,cr,cxb : boolean;
69: jj,ju,kp1 : integer;
70: t,temp : Complex_Number;
71:
72: begin
73: info := 0; -- set info flag
74: cqy := (job/10000 /= 0); -- determine what to compute
75: cqty := (job mod 10000 /= 0);
76: cb := ((job mod 1000)/100 /= 0);
77: cr := ((job mod 100)/10 /= 0);
78: cxb := ((job mod 10) /= 0);
79: ju := min0(k,n-1);
80: if ju = 0 -- special action when n = 1
81: then if cqy then qy(1) := y(1); end if;
82: if cqty then qty(1) := y(1); end if;
83: if cxb then xb(1) := y(1); end if;
84: if cb
85: then if AbsVal(x(1,1)) = 0.0
86: then info := 1;
87: else b(1) := y(1)/x(1,1);
88: end if;
89: end if;
90: if cr then rsd(1) := Create(0.0); end if;
91: return;
92: end if;
93: if cqy -- set up to compute qy or qty
94: then zcopy(n,y'first,y,qy);
95: end if;
96: if cqty
97: then zcopy(n,y'first,y,qty);
98: end if;
99: if cqy -- compute qy
100: then for j in 1..ju loop
101: jj := ju - j + 1;
102: if AbsVal(qraux(jj)) /= 0.0
103: then temp := x(jj,jj);
104: x(jj,jj) := qraux(jj);
105: t := -zdotc(jj,x,qy)/x(jj,jj);
106: zaxpy(n-jj+1,jj,jj,t,x,qy);
107: x(jj,jj) := temp;
108: end if;
109: end loop;
110: end if;
111: if cqty -- compute trans(q)*y
112: then for j in 1..ju loop
113: if AbsVal(qraux(j)) /= 0.0
114: then temp := x(j,j);
115: x(j,j) := qraux(j);
116: t := -zdotc(j,x,qty)/x(j,j);
117: zaxpy(n-j+1,j,j,t,x,qty);
118: x(j,j) := temp;
119: end if;
120: end loop;
121: end if;
122: if cb -- set up to compute b,rsd, or xb
123: then zcopy(k,qty'first,qty,b);
124: end if;
125: kp1 := k + 1;
126: if cxb then zcopy(k,qty'first,qty,xb); end if;
127: if (cr and (k < n))
128: then zcopy(n-k,kp1,qty,rsd);
129: end if;
130: if (cxb and (kp1 <= n))
131: then for i in kp1..n loop
132: xb(i) := Create(0.0);
133: end loop;
134: end if;
135: if cr
136: then for i in 1..k loop
137: rsd(i) := Create(0.0);
138: end loop;
139: end if;
140: if cb -- compute b
141: then for j in 1..k loop
142: jj := k - j + 1;
143: if AbsVal(x(jj,jj)) = 0.0
144: then info := jj; exit;
145: end if;
146: b(jj) := b(jj)/x(jj,jj);
147: if jj /= 1
148: then t := -b(jj);
149: zaxpy(jj-1,1,jj,t,x,b);
150: end if;
151: end loop;
152: end if;
153: if cr or cxb -- compute rsd or xb as requested
154: then for j in 1..ju loop
155: jj := ju - j + 1;
156: if AbsVal(qraux(jj)) /= 0.0
157: then temp := x(jj,jj);
158: x(jj,jj) := qraux(jj);
159: if cr
160: then t := -zdotc(jj,x,rsd)/x(jj,jj);
161: zaxpy(n-jj+1,jj,jj,t,x,rsd);
162: end if;
163: if cxb
164: then t := -zdotc(jj,x,xb)/x(jj,jj);
165: zaxpy(n-jj+1,jj,jj,t,x,xb);
166: end if;
167: x(jj,jj) := temp;
168: end if;
169: end loop;
170: end if;
171: end QRLS;
172:
173: end Standard_Complex_Least_Squares;