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