![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to standard_complex_linear_solvers.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_linear_solvers.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices |
1.1 maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
2:
3: package body Standard_Complex_Linear_Solvers is
4:
5: -- AUXLILIARIES :
6:
7: function cabs ( c : Complex_Number ) return double_float is
8: begin
9: return (ABS(REAL_PART(c)) + ABS(IMAG_PART(c)));
10: end cabs;
11:
12: function dconjg ( x : Complex_Number ) return Complex_Number is
13: begin
14: return Create(REAL_PART(x),-IMAG_PART(x));
15: end dconjg;
16:
17: function csign ( x,y : Complex_Number ) return Complex_Number is
18: begin
19: return (Create(cabs(x)) * y / Create(cabs(y)));
20: end csign;
21:
22: -- TARGET ROUTINES :
23:
24: procedure Scale ( a : in out Matrix; b : in out Vector ) is
25:
26: fac : Complex_Number;
27:
28: function Maximum ( a : in Matrix; i : in integer ) return Complex_Number is
29:
30: res : integer := a'first(2);
31: max : double_float := cabs(a(i,res));
32: tmp : double_float;
33:
34: begin
35: for j in a'first(2)+1..a'last(2) loop
36: tmp := cabs(a(i,j));
37: if tmp > max
38: then max := tmp; res := j;
39: end if;
40: end loop;
41: return a(i,res);
42: end Maximum;
43:
44: procedure Divide ( a : in out Matrix; b : in out Vector;
45: i : in integer; fac : in Complex_Number ) is
46: begin
47: for j in a'range(2) loop
48: a(i,j) := a(i,j)/fac;
49: end loop;
50: b(i) := b(i)/fac;
51: end Divide;
52:
53: begin
54: for i in a'range(1) loop
55: fac := Maximum(a,i);
56: Divide(a,b,i,fac);
57: end loop;
58: end Scale;
59:
60: -- TARGET ROUTINES :
61:
62: procedure lufac ( a : in out Matrix; n : in integer;
63: ipvt : out Standard_Natural_Vectors.Vector;
64: info : out natural ) is
65:
66: kp1,l,nm1 : integer;
67: smax : double_float;
68: temp : Complex_Number;
69:
70: begin
71: info := 0;
72: nm1 := n - 1;
73: if nm1 >= 1
74: then for k in 1..nm1 loop
75: kp1 := k + 1;
76:
77: -- find the pivot index l
78:
79: l := k; smax := cabs(a(k,k)); --modulus(a(k,k));
80: for i in kp1..n loop
81: if cabs(a(i,k)) > smax --modulus(a(i,k)) > smax
82: then l := i;
83: smax := cabs(a(i,k)); --modulus(a(i,k));
84: end if;
85: end loop;
86: ipvt(k) := l;
87:
88: if smax = 0.0
89: then -- this column is already triangularized
90: info := k;
91: else
92:
93: -- interchange if necessary
94:
95: if l /= k
96: then temp := a(l,k);
97: a(l,k) := a(k,k);
98: a(k,k) := temp;
99: end if;
100:
101: -- compute multipliers
102:
103: temp := -Create(1.0)/a(k,k);
104: for i in kp1..n loop
105: a(i,k) := temp * a(i,k);
106: end loop;
107:
108: -- row elimination with column indexing
109:
110: for j in kp1..n loop
111: temp := a(l,j);
112: if l /= k
113: then a(l,j) := a(k,j);
114: a(k,j) := temp;
115: end if;
116: for i in kp1..n loop
117: a(i,j) := a(i,j) + temp * a(i,k);
118: end loop;
119: end loop;
120:
121: end if;
122: end loop;
123: end if;
124: ipvt(n) := n;
125: if AbsVal(a(n,n)) = 0.0
126: then info := n;
127: end if;
128: end lufac;
129:
130: procedure lufco ( a : in out Matrix; n : in integer;
131: ipvt : out Standard_Natural_Vectors.Vector;
132: rcond : out double_float ) is
133:
134: -- NOTE :
135: -- rcond = 1/(norm(a)*(estimate of norm(inverse(a))))
136: -- estimate = norm(z)/norm(y) where a*z = y and ctrans(a)*y = e.
137: -- ctrans(a) is the conjugate transpose of a.
138: -- The components of e are chosen to cause maximum local
139: -- growth in teh elements of w where ctrans(u)*w = e.
140: -- The vectors are frequently rescaled to avoid overflow.
141:
142: z : Standard_Complex_Vectors.Vector(1..n);
143: info,kb,kp1,l : integer;
144: s,sm,sum,anorm,ynorm : double_float;
145: ek,t,wk,wkm : Complex_Number;
146: ipvtt : Standard_Natural_Vectors.Vector(1..n);
147:
148: begin
149: anorm := 0.0; -- compute 1-norm of a
150: for j in 1..n loop
151: sum := 0.0;
152: for i in 1..n loop
153: sum := sum + cabs(a(i,j));
154: end loop;
155: if sum > anorm
156: then anorm := sum;
157: end if;
158: end loop;
159: lufac(a,n,ipvtt,info); -- factor
160: for i in 1..n loop
161: ipvt(i) := ipvtt(i);
162: end loop;
163: ek := Create(1.0); -- solve ctrans(u)*w = e
164: for j in 1..n loop
165: z(j) := Create(0.0);
166: end loop;
167: for k in 1..n loop
168: if cabs(z(k)) /= 0.0
169: then ek := csign(ek,-z(k));
170: end if;
171: if cabs(ek-z(k)) > cabs(a(k,k))
172: then s := cabs(a(k,k))/cabs(ek-z(k));
173: z := Create(s) * z;
174: ek := Create(s) * ek;
175: end if;
176: wk := ek - z(k);
177: wkm := -ek - z(k);
178: s := cabs(wk);
179: sm := cabs(wkm);
180: if cabs(a(k,k)) = 0.0
181: then wk := Create(1.0);
182: wkm := Create(1.0);
183: else wk := wk / dconjg(a(k,k));
184: wkm := wkm / dconjg(a(k,k));
185: end if;
186: kp1 := k + 1;
187: if kp1 <= n
188: then for j in kp1..n loop
189: sm := sm + cabs(z(j)+wkm*dconjg(a(k,j)));
190: z(j) := z(j) + wk*dconjg(a(k,j));
191: s := s + cabs(z(j));
192: end loop;
193: if s < sm
194: then t := wkm - wk;
195: wk := wkm;
196: for j in kp1..n loop
197: z(j) := z(j) + t*dconjg(a(k,j));
198: end loop;
199: end if;
200: end if;
201: z(k) := wk;
202: end loop;
203: sum := 0.0;
204: for i in 1..n loop
205: sum := sum + cabs(z(i));
206: end loop;
207: s := 1.0 / sum;
208: z := Create(s) * z;
209: for k in 1..n loop -- solve ctrans(l)*y = w
210: kb := n+1-k;
211: if kb < n
212: then t := Create(0.0);
213: for i in (kb+1)..n loop
214: t := t + dconjg(a(i,kb))*z(i);
215: end loop;
216: z(kb) := z(kb) + t;
217: end if;
218: if cabs(z(kb)) > 1.0
219: then s := 1.0 / cabs(z(kb));
220: z := Create(s) * z;
221: end if;
222: l := ipvtt(kb);
223: t := z(l);
224: z(l) := z(kb);
225: z(kb) := t;
226: end loop;
227: sum := 0.0;
228: for i in 1..n loop
229: sum := sum + cabs(z(i));
230: end loop;
231: s := 1.0 / sum;
232: z := Create(s) * z;
233: ynorm := 1.0;
234: for k in 1..n loop -- solve l*v = y
235: l := ipvtt(k);
236: t := z(l);
237: z(l) := z(k);
238: z(k) := t;
239: if k < n
240: then for i in (k+1)..n loop
241: z(i) := z(i) + t * a(i,k);
242: end loop;
243: end if;
244: if cabs(z(k)) > 1.0
245: then s := 1.0 / cabs(z(k));
246: z := Create(s) * z;
247: ynorm := s * ynorm;
248: end if;
249: end loop;
250: sum := 0.0;
251: for i in 1..n loop
252: sum := sum + cabs(z(i));
253: end loop;
254: s := 1.0 / sum;
255: z := Create(s) * z;
256: ynorm := s * ynorm;
257: for k in 1..n loop -- solve u*z = v
258: kb := n+1-k;
259: if cabs(z(kb)) > cabs(a(kb,kb))
260: then s := cabs(a(kb,kb)) / cabs(z(kb));
261: z := Create(s) * z;
262: ynorm := s * ynorm;
263: end if;
264: if cabs(a(kb,kb)) = 0.0
265: then z(kb) := Create(1.0);
266: else z(kb) := z(kb) / a(kb,kb);
267: end if;
268: t := -z(kb);
269: for i in 1..(kb-1) loop
270: z(i) := z(i) + t * a(i,kb);
271: end loop;
272: end loop;
273: sum := 0.0; -- make znorm = 1.0
274: for i in 1..n loop
275: sum := sum + cabs(z(i));
276: end loop;
277: s := 1.0 / sum;
278: z := Create(s) * z;
279: ynorm := s * ynorm;
280: if anorm = 0.0
281: then rcond := 0.0;
282: else rcond := ynorm/anorm;
283: end if;
284: end lufco;
285:
286: procedure lusolve ( a : in Matrix; n : in integer;
287: ipvt : in Standard_Natural_Vectors.Vector;
288: b : in out Vector ) is
289:
290: l,nm1,kb : integer;
291: temp : Complex_Number;
292:
293: begin
294: nm1 := n-1;
295: if nm1 >= 1 -- solve l*y = b
296: then for k in 1..nm1 loop
297: l := ipvt(k);
298: temp := b(l);
299: if l /= k
300: then b(l) := b(k);
301: b(k) := temp;
302: end if;
303: for i in (k+1)..n loop
304: b(i) := b(i) + temp * a(i,k);
305: end loop;
306: end loop;
307: end if;
308: for k in 1..n loop -- solve u*x = y
309: kb := n+1-k;
310: b(kb) := b(kb) / a(kb,kb);
311: temp := -b(kb);
312: for j in 1..(kb-1) loop
313: b(j) := b(j) + temp * a(j,kb);
314: end loop;
315: end loop;
316: end lusolve;
317:
318: procedure Triangulate ( a : in out Matrix; n,m : in integer ) is
319:
320: max,cbs : double_float;
321: temp : Complex_Number;
322: pivot,k,kcolumn : integer;
323: tol : constant double_float := 10.0**(-10);
324:
325: begin
326: k := 1;
327: kcolumn := 1;
328: while (k <= n) and (kcolumn <= m) loop
329: max := 0.0; -- find pivot
330: pivot := 0;
331: for l in k..n loop
332: cbs := cabs(a(l,kcolumn));
333: if (cbs > tol) and then (cbs > max)
334: then max := cbs;
335: pivot := l;
336: end if;
337: end loop;
338: if pivot = 0
339: then kcolumn := kcolumn + 1;
340: else if pivot /= k -- interchange if necessary
341: then for i in 1..m loop
342: temp := a(pivot,i);
343: a(pivot,i) := a(k,i);
344: a(k,i) := temp;
345: end loop;
346: end if;
347: for j in (kcolumn+1)..m loop -- triangulate a
348: a(k,j) := a(k,j) / a(k,kcolumn);
349: end loop;
350: a(k,kcolumn) := Create(1.0);
351: for i in (k+1)..n loop
352: for j in (kcolumn+1)..m loop
353: a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
354: end loop;
355: a(i,kcolumn) := Create(0.0);
356: end loop;
357: k := k + 1;
358: kcolumn := kcolumn + 1;
359: end if;
360: end loop;
361: end Triangulate;
362:
363: procedure Diagonalize ( a : in out Matrix; n,m : in integer ) is
364:
365: max : double_float;
366: temp : Complex_Number;
367: pivot,k,kcolumn : integer;
368:
369: begin
370: k := 1;
371: kcolumn := 1;
372: while (k <= n) and (kcolumn <= m) loop
373: max := 0.0; -- find pivot
374: for l in k..n loop
375: if cabs(a(l,kcolumn)) > max
376: then max := cabs(a(l,kcolumn));
377: pivot := l;
378: end if;
379: end loop;
380: if max = 0.0
381: then kcolumn := kcolumn + 1;
382: else if pivot /= k -- interchange if necessary
383: then for i in 1..m loop
384: temp := a(pivot,i);
385: a(pivot,i) := a(k,i);
386: a(k,i) := temp;
387: end loop;
388: end if;
389: for j in (kcolumn+1)..m loop -- diagonalize a
390: a(k,j) := a(k,j) / a(k,kcolumn);
391: end loop;
392: a(k,kcolumn) := Create(1.0);
393: for i in 1..(k-1) loop
394: for j in (kcolumn+1)..m loop
395: a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
396: end loop;
397: end loop;
398: for i in (k+1)..n loop
399: for j in (kcolumn+1)..m loop
400: a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
401: end loop;
402: end loop;
403: for j in 1..(k-1) loop
404: a(j,kcolumn) := Create(0.0);
405: end loop;
406: for j in (k+1)..n loop
407: a(j,kcolumn) := Create(0.0);
408: end loop;
409: k := k + 1;
410: kcolumn := kcolumn + 1;
411: end if;
412: end loop;
413: end Diagonalize;
414:
415: end Standard_Complex_Linear_Solvers;