![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to multprec_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 multprec_complex_linear_solvers.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Matrices |
1.1 maekawa 1: with Multprec_Complex_Numbers; use Multprec_Complex_Numbers;
2:
3: package body Multprec_Complex_Linear_Solvers is
4:
5: -- AUXLILIARIES :
6:
7: function dconjg ( x : Complex_Number ) return Complex_Number is
8:
9: -- DESCRIPTION :
10: -- Returns the complex conjugate of x.
11:
12: res : Complex_Number;
13: re : Floating_Number := REAL_PART(x);
14: im : Floating_Number := IMAG_PART(x);
15:
16: begin
17: Min(im);
18: res := Create(re,im);
19: Clear(re); Clear(im);
20: return res;
21: end dconjg;
22:
23: function csign ( x,y : Complex_Number ) return Complex_Number is
24:
25: -- DESCRIPTION :
26: -- Returns |x|*y/|y|.
27:
28: res : Complex_Number;
29: fltacc : Floating_Number;
30: cmpacc : Complex_Number;
31:
32: begin
33: fltacc := AbsVal(x);
34: res := Create(fltacc);
35: Clear(fltacc);
36: Mul(res,y);
37: fltacc := AbsVal(y);
38: cmpacc := Create(fltacc);
39: Div(res,cmpacc);
40: Clear(fltacc);
41: Clear(cmpacc);
42: return res;
43: end csign;
44:
45: function Inverse_Abs_Sum ( z : Vector ) return Floating_Number is
46:
47: -- DESCRIPTION :
48: -- Returns the reciprocal of the sum of the absolute values in z.
49:
50: res,sum,acc : Floating_Number;
51:
52: begin
53: sum := Create(0);
54: for i in z'range loop
55: acc := AbsVal(z(i));
56: Add(sum,acc);
57: Clear(acc);
58: end loop;
59: res := Create(1);
60: Div(res,sum);
61: Clear(sum);
62: return res;
63: end Inverse_Abs_Sum;
64:
65: -- TARGET ROUTINES :
66:
67: procedure Scale ( a : in out Matrix; b : in out Vector ) is
68:
69: fac : Complex_Number;
70:
71: function Maximum ( a : in Matrix; i : in integer ) return Complex_Number is
72:
73: res : integer := a'first(2);
74: max : Floating_Number := AbsVal(a(i,res));
75: tmp : Floating_Number;
76:
77: begin
78: for j in a'first(2)+1..a'last(2) loop
79: tmp := AbsVal(a(i,j));
80: if tmp > max
81: then max := tmp; res := j;
82: end if;
83: end loop;
84: return a(i,res);
85: end Maximum;
86:
87: procedure Divide ( a : in out Matrix; b : in out Vector;
88: i : in integer; fac : in Complex_Number ) is
89: begin
90: for j in a'range(2) loop
91: a(i,j) := a(i,j)/fac;
92: end loop;
93: b(i) := b(i)/fac;
94: end Divide;
95:
96: begin
97: for i in a'range(1) loop
98: fac := Maximum(a,i);
99: Divide(a,b,i,fac);
100: end loop;
101: end Scale;
102:
103: -- TARGET ROUTINES :
104:
105: procedure lufac ( a : in out Matrix; n : in integer;
106: ipvt : out Standard_Natural_Vectors.Vector;
107: info : out natural ) is
108:
109: kp1,l,nm1 : integer;
110: smax,fltacc : Floating_Number;
111: temp,cmpacc : Complex_Number;
112:
113: begin
114: info := 0;
115: nm1 := n - 1;
116: if nm1 >= 1
117: then for k in 1..nm1 loop
118: kp1 := k + 1;
119: l := k;
120: smax := AbsVal(a(k,k)); -- find the pivot index l
121: for i in kp1..n loop
122: fltacc := AbsVal(a(i,k));
123: if fltacc > smax
124: then l := i;
125: Copy(fltacc,smax);
126: end if;
127: Clear(fltacc);
128: end loop;
129: ipvt(k) := l;
130: if Equal(smax,0.0) -- this column is already triangularized
131: then info := k;
132: else if l /= k -- interchange if necessary
133: then Copy(a(l,k),temp);
134: Copy(a(k,k),a(l,k));
135: Copy(temp,a(k,k)); Clear(temp);
136: end if;
137: cmpacc := Create(1); -- compute multipliers
138: Div(cmpacc,a(k,k));
139: Min(cmpacc);
140: for i in kp1..n loop
141: Mul(a(i,k),cmpacc);
142: end loop;
143: Clear(cmpacc);
144: for j in kp1..n loop -- row elimination
145: Copy(a(l,j),temp);
146: if l /= k
147: then Copy(a(k,j),a(l,j));
148: Copy(temp,a(k,j));
149: end if;
150: for i in kp1..n loop
151: cmpacc := temp*a(i,k);
152: Add(a(i,j),cmpacc);
153: Clear(cmpacc);
154: end loop;
155: end loop;
156: Clear(temp);
157: end if;
158: Clear(smax);
159: end loop;
160: end if;
161: ipvt(n) := n;
162: fltacc := AbsVal(a(n,n));
163: if Equal(fltacc,0.0)
164: then info := n;
165: end if;
166: Clear(fltacc);
167: end lufac;
168:
169: procedure lufco ( a : in out Matrix; n : in integer;
170: ipvt : out Standard_Natural_Vectors.Vector;
171: rcond : out Floating_Number ) is
172:
173: -- NOTE :
174: -- rcond = 1/(norm(a)*(estimate of norm(inverse(a))))
175: -- estimate = norm(z)/norm(y) where a*z = y and ctrans(a)*y = e.
176: -- ctrans(a) is the conjugate transpose of a.
177: -- The components of e are chosen to cause maximum local
178: -- growth in teh elements of w where ctrans(u)*w = e.
179: -- The vectors are frequently rescaled to avoid overflow.
180:
181: z : Multprec_Complex_Vectors.Vector(1..n);
182: info,kb,kp1,l : integer;
183: s,sm,sum,anorm,ynorm,fltacc1,fltacc2 : Floating_Number;
184: ek,t,wk,wkm,cmpacc1,cmpacc2 : Complex_Number;
185: ipvtt : Standard_Natural_Vectors.Vector(1..n);
186:
187: begin
188: anorm := Create(0); -- compute 1-norm of a
189: for j in 1..n loop
190: sum := Create(0);
191: for i in 1..n loop
192: fltacc1 := AbsVal(a(i,j));
193: Add(sum,fltacc1);
194: Clear(fltacc1);
195: end loop;
196: if sum > anorm
197: then Copy(sum,anorm);
198: end if;
199: Clear(sum);
200: end loop;
201: lufac(a,n,ipvtt,info); -- factor
202: ipvt := ipvtt;
203: ek := Create(1); -- solve ctrans(u)*w = e
204: for j in 1..n loop
205: z(j) := Create(0);
206: end loop;
207: for k in 1..n loop
208: fltacc1 := AbsVal(z(k));
209: if not Equal(fltacc1,0.0)
210: then cmpacc1 := -z(k);
211: Clear(ek);
212: ek := csign(ek,cmpacc1);
213: Clear(cmpacc1);
214: end if;
215: Clear(fltacc1);
216: cmpacc1 := ek - z(k);
217: fltacc1 := AbsVal(cmpacc1);
218: Clear(cmpacc1);
219: fltacc2 := AbsVal(a(k,k));
220: if fltacc1 > fltacc2
221: then s := fltacc2/fltacc1;
222: cmpacc1 := Create(s);
223: Mul(z,cmpacc1);
224: Mul(ek,cmpacc1);
225: Clear(cmpacc1);
226: Clear(s);
227: end if;
228: Clear(fltacc1); Clear(fltacc2);
229: Clear(wk); wk := ek - z(k);
230: Clear(wkm); wkm := ek + z(k);
231: Min(wkm);
232: s := AbsVal(wk);
233: sm := AbsVal(wkm);
234: fltacc1 := AbsVal(a(k,k));
235: if Equal(fltacc1,0.0)
236: then Clear(wk); wk := Create(1);
237: Clear(wkm); wkm := Create(1);
238: else cmpacc1 := dconjg(a(k,k));
239: Div(wk,cmpacc1);
240: Div(wkm,cmpacc1);
241: Clear(cmpacc1);
242: end if;
243: Clear(fltacc1);
244: kp1 := k + 1;
245: if kp1 <= n
246: then for j in kp1..n loop
247: cmpacc2 := dconjg(a(k,j));
248: cmpacc1 := wkm*cmpacc2;
249: Add(cmpacc1,z(j));
250: fltacc1 := AbsVal(cmpacc1);
251: Add(sm,fltacc1);
252: Clear(fltacc1);
253: Clear(cmpacc1);
254: cmpacc1 := wk*cmpacc2;
255: Add(z(j),cmpacc1);
256: Clear(cmpacc1);
257: Clear(cmpacc2);
258: fltacc1 := AbsVal(z(j));
259: Add(s,fltacc1);
260: Clear(fltacc1);
261: end loop;
262: if s < sm
263: then t := wkm - wk;
264: Copy(wkm,wk);
265: for j in kp1..n loop
266: cmpacc1 := t*dconjg(a(k,j));
267: Add(z(j),cmpacc1);
268: Clear(cmpacc1);
269: end loop;
270: Clear(t);
271: end if;
272: end if;
273: Copy(wk,z(k));
274: end loop;
275: s := Inverse_Abs_Sum(z);
276: cmpacc1 := Create(s);
277: Mul(z,cmpacc1);
278: Clear(cmpacc1);
279: Clear(s);
280: for k in 1..n loop -- solve ctrans(l)*y = w
281: kb := n+1-k;
282: if kb < n
283: then t := Create(0);
284: for i in (kb+1)..n loop
285: cmpacc1 := dconjg(a(i,kb))*z(i);
286: Add(t,cmpacc1);
287: Clear(cmpacc1);
288: end loop;
289: Add(z(kb),t);
290: end if;
291: fltacc1 := AbsVal(z(kb));
292: if fltacc1 > 1.0
293: then s := Create(1);
294: Div(s,fltacc1);
295: cmpacc1 := Create(s);
296: Mul(z,cmpacc1);
297: Clear(cmpacc1);
298: Clear(s);
299: end if;
300: Clear(fltacc1);
301: l := ipvtt(kb);
302: if l /= kb
303: then Copy(z(l),t);
304: Copy(z(kb),z(l));
305: Copy(t,z(kb));
306: Clear(t);
307: end if;
308: end loop;
309: s := Inverse_Abs_Sum(z);
310: cmpacc1 := Create(s);
311: Clear(s);
312: Mul(z,cmpacc1);
313: Clear(cmpacc1);
314: ynorm := Create(1);
315: for k in 1..n loop -- solve l*v = y
316: l := ipvtt(k);
317: if l /= k
318: then Copy(z(l),t);
319: Copy(z(k),z(l));
320: Copy(t,z(k));
321: else Copy(z(l),t);
322: end if;
323: if k < n
324: then for i in (k+1)..n loop
325: cmpacc1 := t*a(i,k);
326: Add(z(i),cmpacc1);
327: Clear(cmpacc1);
328: end loop;
329: end if;
330: Clear(t);
331: fltacc1 := AbsVal(z(k));
332: if fltacc1 > 1.0
333: then s := Create(1);
334: Div(s,fltacc1);
335: cmpacc1 := Create(s);
336: Mul(z,cmpacc1);
337: Clear(cmpacc1);
338: Mul(ynorm,s);
339: Clear(s);
340: end if;
341: Clear(fltacc1);
342: end loop;
343: s := Inverse_Abs_Sum(z);
344: cmpacc1 := Create(s);
345: Mul(z,cmpacc1);
346: Clear(cmpacc1);
347: Mul(ynorm,s);
348: Clear(s);
349: for k in 1..n loop -- solve u*z = v
350: kb := n+1-k;
351: fltacc1 := AbsVal(z(kb));
352: fltacc2 := AbsVal(a(kb,kb));
353: if fltacc1 > fltacc2
354: then s := fltacc2/fltacc1;
355: cmpacc1 := Create(s);
356: Mul(z,cmpacc1);
357: Mul(ynorm,s);
358: Clear(s);
359: end if;
360: if Equal(fltacc2,0.0)
361: then Clear(z(kb));
362: z(kb) := Create(1);
363: else Div(z(kb),a(kb,kb));
364: end if;
365: Clear(fltacc1);
366: Clear(fltacc2);
367: t := -z(kb);
368: for i in 1..(kb-1) loop
369: cmpacc1 := t*a(i,kb);
370: Add(z(i),cmpacc1);
371: Clear(cmpacc1);
372: end loop;
373: Clear(t);
374: end loop;
375: s := Inverse_Abs_Sum(z); -- make znorm = 1.0
376: cmpacc1 := Create(s);
377: Mul(z,cmpacc1);
378: Clear(cmpacc1);
379: Mul(ynorm,s);
380: if Equal(anorm,0.0)
381: then rcond := Create(0);
382: else rcond := ynorm/anorm;
383: end if;
384: Clear(ynorm); Clear(s); Clear(sm);
385: Clear(ek); Clear(t); Clear(wk); Clear(wkm);
386: Clear(z);
387: end lufco;
388:
389: procedure lusolve ( a : in Matrix; n : in integer;
390: ipvt : in Standard_Natural_Vectors.Vector;
391: b : in out Vector ) is
392:
393: l,nm1,kb : integer;
394: temp,acc : Complex_Number;
395:
396: begin
397: nm1 := n-1;
398: if nm1 >= 1 -- solve l*y = b
399: then for k in 1..nm1 loop
400: l := ipvt(k);
401: Copy(b(l),temp);
402: if l /= k
403: then Copy(b(k),b(l));
404: Copy(temp,b(k));
405: end if;
406: for i in (k+1)..n loop
407: acc := temp*a(i,k);
408: Add(b(i),acc);
409: Clear(acc);
410: end loop;
411: Clear(temp);
412: end loop;
413: end if;
414: for k in 1..n loop -- solve u*x = y
415: kb := n+1-k;
416: Div(b(kb),a(kb,kb));
417: temp := -b(kb);
418: for j in 1..(kb-1) loop
419: acc := temp*a(j,kb);
420: Add(b(j),acc);
421: Clear(acc);
422: end loop;
423: Clear(temp);
424: end loop;
425: end lusolve;
426:
427: procedure Triangulate ( a : in out Matrix; n,m : in integer ) is
428:
429: max,cbs : Floating_Number;
430: temp : Complex_Number;
431: pivot,k,kcolumn : integer;
432:
433: begin
434: k := 1;
435: kcolumn := 1;
436: while (k <= n) and (kcolumn <= m) loop
437: max := Create(0); -- find pivot
438: pivot := 0;
439: for l in k..n loop
440: cbs := AbsVal(a(l,kcolumn));
441: if cbs > max
442: then max := cbs;
443: pivot := l;
444: end if;
445: end loop;
446: if pivot = 0
447: then kcolumn := kcolumn + 1;
448: else if pivot /= k -- interchange if necessary
449: then for i in 1..m loop
450: temp := a(pivot,i);
451: a(pivot,i) := a(k,i);
452: a(k,i) := temp;
453: end loop;
454: end if;
455: for j in (kcolumn+1)..m loop -- triangulate a
456: a(k,j) := a(k,j) / a(k,kcolumn);
457: end loop;
458: a(k,kcolumn) := Create(1);
459: for i in (k+1)..n loop
460: for j in (kcolumn+1)..m loop
461: a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
462: end loop;
463: a(i,kcolumn) := Create(0);
464: end loop;
465: k := k + 1;
466: kcolumn := kcolumn + 1;
467: end if;
468: end loop;
469: end Triangulate;
470:
471: procedure Diagonalize ( a : in out Matrix; n,m : in integer ) is
472:
473: max : Floating_Number;
474: temp : Complex_Number;
475: pivot,k,kcolumn : integer;
476:
477: begin
478: k := 1;
479: kcolumn := 1;
480: while (k <= n) and (kcolumn <= m) loop
481: max := Create(0); -- find pivot
482: for l in k..n loop
483: if AbsVal(a(l,kcolumn)) > max
484: then max := AbsVal(a(l,kcolumn));
485: pivot := l;
486: end if;
487: end loop;
488: if Equal(max,0.0)
489: then kcolumn := kcolumn + 1;
490: else if pivot /= k -- interchange if necessary
491: then for i in 1..m loop
492: temp := a(pivot,i);
493: a(pivot,i) := a(k,i);
494: a(k,i) := temp;
495: end loop;
496: end if;
497: for j in (kcolumn+1)..m loop -- diagonalize a
498: a(k,j) := a(k,j) / a(k,kcolumn);
499: end loop;
500: a(k,kcolumn) := Create(1);
501: for i in 1..(k-1) loop
502: for j in (kcolumn+1)..m loop
503: a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
504: end loop;
505: end loop;
506: for i in (k+1)..n loop
507: for j in (kcolumn+1)..m loop
508: a(i,j) := a(i,j) - a(i,kcolumn) * a(k,j);
509: end loop;
510: end loop;
511: for j in 1..(k-1) loop
512: a(j,kcolumn) := Create(0);
513: end loop;
514: for j in (k+1)..n loop
515: a(j,kcolumn) := Create(0);
516: end loop;
517: k := k + 1;
518: kcolumn := kcolumn + 1;
519: end if;
520: end loop;
521: end Diagonalize;
522:
523: end Multprec_Complex_Linear_Solvers;