![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to standard_floating_qr_decomposition.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_qr_decomposition.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_Mathematical_Functions; use Standard_Mathematical_Functions;
3:
4: package body Standard_Floating_QR_Decomposition 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: function dmax1 ( a,b : double_float ) return double_float is
20:
21: -- DESCRIPTION : returns the maximum of a and b.
22:
23: begin
24: if a >= b
25: then return a;
26: else return b;
27: end if;
28: end dmax1;
29:
30: function dsign ( a,b : double_float ) return double_float is
31:
32: -- DESCRIPTION : returns the absolute value of a, set to the same sign as b.
33:
34: begin
35: if b >= 0.0
36: then return abs(a);
37: else return -abs(a);
38: end if;
39: end dsign;
40:
41: procedure dswap ( a : in out Standard_Floating_Matrices.Matrix;
42: k1,k2 : in natural ) is
43:
44: -- DESCRIPTION :
45: -- Swaps the columns k1 and k2 in the matrix a.
46:
47: tmp : double_float;
48:
49: begin
50: for i in a'range(1) loop
51: tmp := a(i,k1); a(i,k1) := a(i,k2); a(i,k2) := tmp;
52: end loop;
53: end dswap;
54:
55: function dnrm2 ( a : Standard_Floating_Matrices.Matrix; row,col : natural )
56: return double_float is
57:
58: -- DESCRIPTION :
59: -- Computes the 2-norm of the vector in the column col of the matrix,
60: -- starting at the given row.
61:
62: sum : double_float := 0.0;
63:
64: begin
65: for i in row..a'last(1) loop
66: sum := sum + a(i,col)*a(i,col);
67: end loop;
68: return SQRT(sum);
69: end dnrm2;
70:
71: function ddot ( a : Standard_Floating_Matrices.Matrix; row,k1,k2 : natural )
72: return double_float is
73:
74: -- DESCRIPTION :
75: -- Returns the inner product of the vectors in the columns k1 and k2,
76: -- starting at the given row.
77:
78: res : double_float := 0.0;
79:
80: begin
81: for i in row..a'last(1) loop
82: res := res + a(i,k1)*a(i,k2);
83: end loop;
84: return res;
85: end ddot;
86:
87: procedure daxpy ( a : in out Standard_Floating_Matrices.Matrix;
88: f : in double_float; row,k1,k2 : in integer ) is
89:
90: -- DESCRIPTION :
91: -- The column k2 is added with f times the column k1, starting at row.
92:
93: begin
94: for i in row..a'last(1) loop
95: a(i,k2) := a(i,k2) + f*a(i,k1);
96: end loop;
97: end daxpy;
98:
99: procedure dscal ( a : in out Standard_Floating_Matrices.Matrix;
100: f : in double_float; row,col : in integer ) is
101:
102: -- DESCRIPTION :
103: -- Multiplies the column col of the matrix with f, starting at row.
104:
105: begin
106: for i in row..a'last(1) loop
107: a(i,col) := f*a(i,col);
108: end loop;
109: end dscal;
110:
111: procedure QRD ( x : in out Standard_Floating_Matrices.Matrix;
112: qraux : in out Standard_Floating_Vectors.Vector;
113: jpvt : in out Standard_Natural_Vectors.Vector;
114: piv : in boolean ) is
115:
116: n : constant natural := x'length(1); -- number of rows
117: p : constant natural := x'length(2); -- number of columns
118: work : Standard_Floating_Vectors.Vector(x'range(2));
119: jj,jp,lp1,lup,maxj,pl,pu : integer;
120: maxnrm,tt,nrmxl,t : double_float;
121: negj,swapj : boolean;
122:
123: begin
124: pl := 1; pu := 0;
125: if piv
126: then for j in x'range(2) loop -- rearrange columns according to jpvt
127: swapj := (jpvt(j) > 0);
128: negj := (jpvt(j) < 0);
129: jpvt(j) := j;
130: if negj
131: then jpvt(j) := -j;
132: end if;
133: if (swapj and then (j /= pl))
134: then dswap(x,pl,j);
135: jpvt(j) := jpvt(pl);
136: jpvt(pl) := j;
137: pl := pl + 1;
138: end if;
139: end loop;
140: pu := p;
141: for j in 1..p loop
142: jj := p - j + 1;
143: if jpvt(jj) < 0
144: then jpvt(jj) := -jpvt(jj);
145: if jj /= pu
146: then dswap(x,pu,jj);
147: jp := jpvt(pu);
148: jpvt(pu) := jpvt(jj);
149: jpvt(jj) := jp;
150: end if;
151: pu := pu - 1;
152: end if;
153: end loop;
154: end if;
155: for j in pl..pu loop -- compute norms of the free columns
156: qraux(j) := dnrm2(x,1,j);
157: work(j) := qraux(j);
158: end loop;
159: lup := min0(n,p); -- perform the householder reduction of x
160: for l in 1..lup loop
161: if (l >= pl and l < pu)
162: then maxnrm := 0.0; -- locate column with largest norm and
163: maxj := l; -- bring it in the pivot position
164: for j in l..pu loop
165: if qraux(j) > maxnrm
166: then maxnrm := qraux(j);
167: maxj := j;
168: end if;
169: end loop;
170: if maxj /= l
171: then dswap(x,l,maxj);
172: qraux(maxj) := qraux(l);
173: work(maxj) := work(l);
174: jp := jpvt(maxj);
175: jpvt(maxj) := jpvt(l);
176: jpvt(l) := jp;
177: end if;
178: end if;
179: qraux(l) := 0.0;
180: if l /= n
181: then nrmxl := dnrm2(x,l,l); -- householder transformation for column l
182: if nrmxl /= 0.0
183: then if x(l,l) /= 0.0
184: then nrmxl := dsign(nrmxl,x(l,l));
185: end if;
186: dscal(x,1.0/nrmxl,l,l);
187: x(l,l) := 1.0 + x(l,l);
188: lp1 := l + 1; -- apply the transformation to the remaining
189: for j in lp1..p loop -- columns, updating the norms
190: t := -ddot(x,l,l,j)/x(l,l);
191: daxpy(x,t,l,l,j);
192: if (j >= pl) and (j <= pu) and (qraux(j) /= 0.0)
193: then tt := 1.0 - (abs(x(l,j))/qraux(j))**2;
194: tt := dmax1(tt,0.0);
195: t := tt;
196: tt := 1.0 + 0.05*tt*(qraux(j)/work(j))**2;
197: if tt /= 1.0
198: then qraux(j) := qraux(j)*SQRT(t);
199: else qraux(j) := dnrm2(x,l+1,j);
200: work(j) := qraux(j);
201: end if;
202: end if;
203: end loop;
204: qraux(l) := x(l,l); -- save the transformation
205: x(l,l) := -nrmxl;
206: end if;
207: end if;
208: end loop;
209: end QRD;
210:
211: procedure Permute_Columns ( x : in out Standard_Floating_Matrices.Matrix;
212: jpvt : in Standard_Natural_Vectors.Vector ) is
213:
214: res : Standard_Floating_Matrices.Matrix(x'range(1),x'range(2));
215:
216: begin
217: for k in jpvt'range loop
218: for i in res'range(1) loop
219: res(i,k) := x(i,jpvt(k));
220: end loop;
221: end loop;
222: x := res;
223: end Permute_Columns;
224:
225: procedure Permute ( x : in out Standard_Floating_Vectors.Vector;
226: jpvt : in Standard_Natural_Vectors.Vector ) is
227:
228: res : Standard_Floating_Vectors.Vector(x'range);
229:
230: begin
231: for k in jpvt'range loop
232: res(k) := x(jpvt(k));
233: end loop;
234: x := res;
235: end Permute;
236:
237: procedure Basis ( qr : in out Standard_Floating_Matrices.Matrix;
238: x : in Standard_Floating_Matrices.Matrix ) is
239:
240: sum : double_float;
241: wrk : Standard_Floating_Vectors.Vector(qr'range(2));
242:
243: begin
244: for j in x'range(2) loop -- compute jth column of q
245: for i in qr'range(1) loop
246: sum := x(i,j);
247: for k in qr'first(2)..(j-1) loop
248: sum := sum - qr(i,k)*qr(k,j);
249: end loop;
250: wrk(i) := sum/qr(j,j);
251: end loop;
252: for i in qr'range(1) loop
253: qr(i,j) := wrk(i);
254: end loop;
255: end loop;
256: end Basis;
257:
258: end Standard_Floating_QR_Decomposition;