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