Annotation of OpenXM_contrib/PHC/Ada/Continuation/correctors.adb, Revision 1.1.1.1
1.1 maekawa 1: with Standard_Natural_Vectors;
2: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
3: with Standard_Complex_Linear_Solvers; use Standard_Complex_Linear_Solvers;
4: with Standard_Complex_QR_Decomposition; use Standard_Complex_QR_Decomposition;
5: with Standard_Complex_Least_Squares; use Standard_Complex_Least_Squares;
6: with Process_io; use Process_io;
7:
8: package body Correctors is
9:
10: -- AUXILIARIES FOR IMPLEMENTING MULTIPLE CORRECTORS :
11:
12: procedure Equals ( s : in Solu_Info_Array; x : in Vector; i : in natural;
13: d : in double_float; j : in out natural ) is
14: eq : boolean;
15:
16: begin
17: while j < i loop
18: eq := true;
19: for k in x'range loop
20: eq := Equal(s(j).sol.v(k),x(k),d);
21: exit when not eq;
22: end loop;
23: exit when eq;
24: j := j + 1;
25: end loop;
26: end Equals;
27:
28: generic
29: with procedure Corrector ( s : in out Solu_Info; c : in Corr_Pars );
30: procedure Multiple_Silent_Corrector
31: ( s : in out Solu_Info_Array;
32: pivot : in out natural; dist_sols : in double_float;
33: c : in Corr_Pars; fail : out boolean );
34:
35: -- DESCRIPTION :
36: -- Allows to build any multiple silent corrector,
37: -- depending on the corrector supplied as generic parameter.
38:
39: generic
40: with procedure Corrector ( file : in file_type;
41: s : in out Solu_Info; c : in Corr_Pars );
42: procedure Multiple_Reporting_Corrector
43: ( file : in file_type; s : in out Solu_Info_Array;
44: pivot : in out natural; dist_sols : in double_float;
45: c : in Corr_Pars; fail : out boolean );
46:
47: -- DESCRIPTION :
48: -- Allows to build any multiple reporting corrector,
49: -- depending on the corrector supplied as generic parameter.
50:
51: procedure Multiple_Silent_Corrector
52: ( s : in out Solu_Info_Array;
53: pivot : in out natural; dist_sols : in double_float;
54: c : in Corr_Pars; fail : out boolean ) is
55:
56: i,j : natural;
57: ffail : boolean := false;
58:
59: begin
60: i := pivot;
61: loop
62: Corrector(s(i),c);
63: ffail := (((s(i).resa > c.epsaf) and then (s(i).cora > c.epsax))
64: or else ((s(i).resr > c.epsrf) and then (s(i).corr > c.epsrx)));
65: if not ffail
66: then j := 1;
67: Equals(s,s(i).sol.v,i,dist_sols,j);
68: if j /= i
69: then ffail := true;
70: end if;
71: end if;
72: exit when ffail;
73: i := i + 1;
74: if i > s'last
75: then i := s'first;
76: end if;
77: exit when (i = pivot);
78: end loop;
79: if ffail
80: then pivot := i;
81: end if;
82: fail := ffail;
83: end Multiple_Silent_Corrector;
84:
85: procedure Multiple_Reporting_Corrector
86: ( file : in file_type; s : in out Solu_Info_Array;
87: pivot : in out natural; dist_sols : in double_float;
88: c : in Corr_Pars; fail : out boolean ) is
89:
90: i,j : natural;
91: ffail : boolean := false;
92:
93: begin
94: i := pivot;
95: loop
96: Write_path(file,i);
97: Corrector(file,s(i),c);
98: sWrite(file,s(i).sol.all);
99: ffail := (((s(i).resa > c.epsaf) and then (s(i).cora > c.epsax))
100: or else ((s(i).resr > c.epsrf) and then (s(i).corr > c.epsrx)));
101: if not ffail
102: then j := 1;
103: Equals(s,s(i).sol.v,i,dist_sols,j);
104: if j /= i
105: then ffail := true;
106: end if;
107: end if;
108: exit when ffail;
109: i := i + 1;
110: if i > s'last
111: then i := s'first;
112: end if;
113: exit when (i = pivot);
114: end loop;
115: if ffail
116: then pivot := i;
117: end if;
118: fail := ffail;
119: end Multiple_Reporting_Corrector;
120:
121: -- TARGET PROCEDURES :
122:
123: procedure Affine_Single_Loose_Normal_Silent_Corrector
124: ( s : in out Solu_Info; c : in Corr_Pars ) is
125:
126: info : integer;
127: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
128: j : Matrix(y'range,y'range);
129: ipvt : Standard_Natural_Vectors.Vector(y'range);
130: nit : natural := 0;
131: normv : double_float;
132:
133: begin
134: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
135: normv := Norm(s.sol.v);
136: if normv + 1.0 /= 1.0
137: then s.corr := s.cora / normv; s.resr := s.resa / normv;
138: end if;
139: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
140: j := dH(s.sol.v,s.sol.t);
141: lufac(j,y'last,ipvt,info);
142: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
143: Min(y);
144: lusolve(j,y'last,ipvt,y); s.cora := Norm(y);
145: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
146: y := H(s.sol.v,s.sol.t); s.resa := Norm(y);
147: normv := Norm(s.sol.v);
148: if normv + 1.0 /= 1.0
149: then s.corr := s.cora / normv; s.resr := s.resa / normv;
150: end if;
151: nit := nit + 1;
152: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
153: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
154: end loop; -- STOP WHEN DESIRED PRECISION REACHED
155: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
156: exception
157: when numeric_error => return;
158: end Affine_Single_Loose_Normal_Silent_Corrector;
159:
160: procedure Affine_Single_Loose_Normal_Reporting_Corrector
161: ( file : in file_type;
162: s : in out Solu_Info; c : in Corr_Pars ) is
163:
164: info : integer;
165: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
166: j : Matrix(y'range,y'range);
167: ipvt : Standard_Natural_Vectors.Vector(y'range);
168: nit : natural := 0;
169: normv : double_float;
170:
171: begin
172: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
173: normv := Norm(s.sol.v);
174: if normv + 1.0 /= 1.0
175: then s.corr := s.cora / normv;
176: s.resr := s.resa / normv;
177: end if;
178: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
179: j := dH(s.sol.v,s.sol.t);
180: lufac(j,y'last,ipvt,info);
181: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
182: Min(y);
183: lusolve(j,y'last,ipvt,y); s.cora := Norm(y);
184: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
185: y := H(s.sol.v,s.sol.t); s.resa := Norm(y);
186: normv := Norm(s.sol.v);
187: if normv + 1.0 /= 1.0
188: then s.corr := s.cora / normv; s.resr := s.resa / normv;
189: end if;
190: cWrite(file,s.cora,s.corr,s.resa,s.resr); -- WRITE STATISTICS
191: nit := nit + 1;
192: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
193: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
194: end loop; -- STOP WHEN DESIRED PRECISION REACHED
195: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
196: exception
197: when numeric_error => return;
198: end Affine_Single_Loose_Normal_Reporting_Corrector;
199:
200: procedure Affine_Single_Loose_Conditioned_Silent_Corrector
201: ( s : in out Solu_Info; c : in Corr_Pars ) is
202:
203: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
204: n : constant natural := y'last;
205: j : Matrix(y'range,y'range);
206: qraux : Standard_Complex_Vectors.Vector(y'range)
207: := (y'range => Create(0.0));
208: dum : Standard_Complex_Vectors.Vector(y'range);
209: ipvt : Standard_Natural_Vectors.Vector(y'range);
210: info : integer;
211: nit : natural := 0;
212: normv : double_float;
213:
214: begin
215: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
216: normv := Norm(s.sol.v);
217: if normv + 1.0 /= 1.0
218: then s.corr := s.cora / normv; s.resr := s.resa / normv;
219: end if;
220: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
221: j := dH(s.sol.v,s.sol.t);
222: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
223: Min(y); -- FOLLOWED BY LEAST SQUARES
224: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
225: s.cora := Norm(y);
226: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
227: y := H(s.sol.v,s.sol.t); s.resa := Norm(y);
228: normv := Norm(s.sol.v);
229: if normv + 1.0 /= 1.0
230: then s.corr := s.cora / normv; s.resr := s.resa / normv;
231: end if;
232: nit := nit + 1;
233: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
234: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
235: end loop; -- STOP WHEN DESIRED PRECISION REACHED
236: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
237: j := dH(s.sol.v,s.sol.t);
238: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
239: exception
240: when numeric_error => return;
241: end Affine_Single_Loose_Conditioned_Silent_Corrector;
242:
243: procedure Affine_Single_Loose_Conditioned_Reporting_Corrector
244: ( file : in file_type;
245: s : in out Solu_Info; c : in Corr_Pars ) is
246:
247: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
248: n : constant natural := y'last;
249: j : Matrix(y'range,y'range);
250: qraux : Standard_Complex_Vectors.Vector(y'range)
251: := (y'range => Create(0.0));
252: dum : Standard_Complex_Vectors.Vector(y'range);
253: ipvt : Standard_Natural_Vectors.Vector(y'range);
254: info : integer;
255: nit : natural := 0;
256: normv : double_float;
257:
258: begin
259: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
260: normv := Norm(s.sol.v);
261: if normv + 1.0 /= 1.0
262: then s.corr := s.cora / normv; s.resr := s.resa / normv;
263: end if;
264: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
265: j := dH(s.sol.v,s.sol.t);
266: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
267: Min(y); -- FOLLOWED BY LEAST SQUARES
268: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
269: s.cora := Norm(y);
270: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
271: y := H(s.sol.v,s.sol.t); s.resa := Norm(y);
272: normv := Norm(s.sol.v);
273: if normv + 1.0 /= 1.0
274: then s.corr := s.cora / normv; s.resr := s.resa / normv;
275: end if;
276: cWrite(file,s.cora,s.corr,s.resa,s.resr); -- WRITE STATISTICS
277: cWrite(file,s.rcond,s.sol.m);
278: nit := nit + 1;
279: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
280: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
281: end loop; -- STOP WHEN DESIRED PRECISION REACHED
282: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
283: j := dH(s.sol.v,s.sol.t);
284: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
285: exception
286: when numeric_error => return;
287: end Affine_Single_Loose_Conditioned_Reporting_Corrector;
288:
289: procedure Affine_Single_Severe_Normal_Silent_Corrector
290: ( s : in out Solu_Info; c : in Corr_Pars ) is
291:
292: info : integer;
293: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
294: j : Matrix(y'range,y'range);
295: ipvt : Standard_Natural_Vectors.Vector(y'range);
296: nit : natural := 0;
297: normv,ncora,nresa,ncorr,nresr : double_float;
298: stop : boolean;
299:
300: begin
301: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
302: normv := Norm(s.sol.v);
303: if normv + 1.0 /= 1.0
304: then s.corr := s.cora / normv; s.resr := s.resa / normv;
305: end if;
306: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
307: j := dH(s.sol.v,s.sol.t);
308: lufac(j,y'last,ipvt,info);
309: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
310: Min(y);
311: lusolve(j,y'last,ipvt,y); ncora := Norm(y);
312: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
313: y := H(s.sol.v,s.sol.t); nresa := Norm(y);
314: normv := Norm(s.sol.v);
315: if normv + 1.0 /= 1.0
316: then ncorr := ncora / normv; nresr := nresa / normv;
317: end if;
318: nit := nit + 1;
319: stop := (((ncora > s.cora) and then (nresa > s.resa))
320: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
321: -- STOP WHEN DIVERGENCE
322: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
323: exit when stop;
324: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
325: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
326: end loop; -- STOP WHEN DESIRED PRECISION REACHED
327: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
328: exception
329: when numeric_error => return;
330: end Affine_Single_Severe_Normal_Silent_Corrector;
331:
332: procedure Affine_Single_Severe_Normal_Reporting_Corrector
333: ( file : in file_type;
334: s : in out Solu_Info; c : in Corr_Pars ) is
335:
336: info : integer;
337: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
338: j : Matrix(y'range,y'range);
339: ipvt : Standard_Natural_Vectors.Vector(y'range);
340: nit : natural := 0;
341: normv,ncora,nresa,ncorr,nresr : double_float;
342: stop : boolean;
343:
344: begin
345: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
346: normv := Norm(s.sol.v);
347: if normv + 1.0 /= 1.0
348: then s.corr := s.cora / normv; s.resr := s.resa / normv;
349: end if;
350: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
351: j := dH(s.sol.v,s.sol.t);
352: lufac(j,y'last,ipvt,info);
353: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
354: Min(y);
355: lusolve(j,y'last,ipvt,y); ncora := Norm(y);
356: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
357: y := H(s.sol.v,s.sol.t); nresa := Norm(y);
358: normv := Norm(s.sol.v);
359: if normv + 1.0 /= 1.0
360: then ncorr := ncora / normv; nresr := nresa / normv;
361: end if;
362: cWrite(file,ncora,ncorr,nresa,nresr); -- WRITE STATISTICS
363: nit := nit + 1;
364: stop := (((ncora > s.cora) and then (nresa > s.resa))
365: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
366: -- STOP WHEN DIVERGENCE
367: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
368: exit when stop;
369: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
370: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
371: end loop; -- STOP WHEN DESIRED PRECISION REACHED
372: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
373: exception
374: when numeric_error => return;
375: end Affine_Single_Severe_Normal_Reporting_Corrector;
376:
377: procedure Affine_Single_Severe_Conditioned_Silent_Corrector
378: ( s : in out Solu_Info; c : in Corr_Pars ) is
379:
380: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
381: n : constant natural := y'last;
382: j : Matrix(y'range,y'range);
383: qraux : Standard_Complex_Vectors.Vector(y'range)
384: := (y'range => Create(0.0));
385: dum : Standard_Complex_Vectors.Vector(y'range);
386: ipvt : Standard_Natural_Vectors.Vector(y'range);
387: info : integer;
388: nit : natural := 0;
389: normv,ncora,nresa,ncorr,nresr : double_float;
390: stop : boolean;
391:
392: begin
393: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
394: normv := Norm(s.sol.v);
395: if normv + 1.0 /= 1.0
396: then s.corr := s.cora / normv; s.resr := s.resa / normv;
397: end if;
398: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
399: j := dH(s.sol.v,s.sol.t);
400: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
401: Min(y); -- FOLLOWED BY LEAST SQUARES
402: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
403: ncora := Norm(y);
404: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
405: y := H(s.sol.v,s.sol.t); nresa := Norm(y);
406: normv := Norm(s.sol.v);
407: if normv + 1.0 /= 1.0
408: then ncorr := ncora / normv; nresr := nresa / normv;
409: end if;
410: nit := nit + 1;
411: stop := (((ncora > s.cora) and then (nresa > s.resa))
412: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
413: -- STOP WHEN DIVERGENCE
414: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
415: exit when stop;
416: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
417: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
418: end loop; -- STOP WHEN DESIRED PRECISION REACHED
419: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
420: j := dH(s.sol.v,s.sol.t);
421: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
422: exception
423: when numeric_error => return;
424: end Affine_Single_Severe_Conditioned_Silent_Corrector;
425:
426: procedure Affine_Single_Severe_Conditioned_Reporting_Corrector
427: ( file : in file_type;
428: s : in out Solu_Info; c : in Corr_Pars ) is
429:
430: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
431: n : constant natural := y'last;
432: j : Matrix(y'range,y'range);
433: qraux : Standard_Complex_Vectors.Vector(y'range)
434: := (y'range => Create(0.0));
435: dum : Standard_Complex_Vectors.Vector(y'range);
436: ipvt : Standard_Natural_Vectors.Vector(y'range);
437: info : integer;
438: nit : natural := 0;
439: normv,ncora,nresa,ncorr,nresr : double_float;
440: stop : boolean;
441:
442: begin
443: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
444: normv := Norm(s.sol.v);
445: if normv + 1.0 /= 1.0
446: then s.corr := s.cora / normv; s.resr := s.resa / normv;
447: end if;
448: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
449: j := dH(s.sol.v,s.sol.t);
450: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
451: Min(y); -- FOLLOWED BY LEAST SQUARES
452: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
453: ncora := Norm(y);
454: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
455: y := H(s.sol.v,s.sol.t); nresa := Norm(y);
456: normv := Norm(s.sol.v);
457: if normv + 1.0 /= 1.0
458: then ncorr := ncora / normv; nresr := nresa / normv;
459: end if;
460: cWrite(file,ncora,ncorr,nresa,nresr); -- WRITE STATISTICS
461: cWrite(file,s.rcond,s.sol.m);
462: nit := nit + 1;
463: stop := (((ncora > s.cora) and then (nresa > s.resa))
464: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
465: -- STOP WHEN DIVERGENCE
466: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
467: exit when stop;
468: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
469: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
470: end loop; -- STOP WHEN DESIRED PRECISION REACHED
471: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
472: j := dH(s.sol.v,s.sol.t);
473: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
474: exception
475: when numeric_error => return;
476: end Affine_Single_Severe_Conditioned_Reporting_Corrector;
477:
478: procedure Affine_Multiple_Loose_Normal_Silent_Corrector
479: ( s : in out Solu_Info_Array;
480: pivot : in out natural; dist_sols : in double_float;
481: c : in Corr_Pars; fail : out boolean ) is
482:
483: procedure Single_Corrector is
484: new Affine_Single_Loose_Normal_Silent_Corrector(Norm,H,dH);
485: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
486:
487: begin
488: Corrector(s,pivot,dist_sols,c,fail);
489: end Affine_Multiple_Loose_Normal_Silent_Corrector;
490:
491: procedure Affine_Multiple_Loose_Normal_Reporting_Corrector
492: ( file : in file_type; s : in out Solu_Info_Array;
493: pivot : in out natural; dist_sols : in double_float;
494: c : in Corr_Pars; fail : out boolean ) is
495:
496: procedure Single_Corrector is
497: new Affine_Single_Loose_Normal_Reporting_Corrector(Norm,H,dH);
498: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
499:
500: begin
501: Corrector(file,s,pivot,dist_sols,c,fail);
502: end Affine_Multiple_Loose_Normal_Reporting_Corrector;
503:
504: procedure Affine_Multiple_Loose_Conditioned_Silent_Corrector
505: ( s : in out Solu_Info_Array;
506: pivot : in out natural; dist_sols : in double_float;
507: c : in Corr_Pars; fail : out boolean ) is
508:
509: procedure Single_Corrector is
510: new Affine_Single_Loose_Conditioned_Silent_Corrector(Norm,H,dH);
511: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
512:
513: begin
514: Corrector(s,pivot,dist_sols,c,fail);
515: end Affine_Multiple_Loose_Conditioned_Silent_Corrector;
516:
517: procedure Affine_Multiple_Loose_Conditioned_Reporting_Corrector
518: ( file : in file_type; s : in out Solu_Info_Array;
519: pivot : in out natural; dist_sols : in double_float;
520: c : in Corr_Pars; fail : out boolean ) is
521:
522: procedure Single_Corrector is
523: new Affine_Single_Loose_Conditioned_Reporting_Corrector(Norm,H,dH);
524: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
525:
526: begin
527: Corrector(file,s,pivot,dist_sols,c,fail);
528: end Affine_Multiple_Loose_Conditioned_Reporting_Corrector;
529:
530: procedure Affine_Multiple_Severe_Normal_Silent_Corrector
531: ( s : in out Solu_Info_Array;
532: pivot : in out natural; dist_sols : in double_float;
533: c : in Corr_Pars; fail : out boolean ) is
534:
535: procedure Single_Corrector is
536: new Affine_Single_Severe_Normal_Silent_Corrector(Norm,H,dH);
537: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
538:
539: begin
540: Corrector(s,pivot,dist_sols,c,fail);
541: end Affine_Multiple_Severe_Normal_Silent_Corrector;
542:
543: procedure Affine_Multiple_Severe_Normal_Reporting_Corrector
544: ( file : in file_type; s : in out Solu_Info_Array;
545: pivot : in out natural; dist_sols : in double_float;
546: c : in Corr_Pars; fail : out boolean ) is
547:
548: procedure Single_Corrector is
549: new Affine_Single_Severe_Normal_Reporting_Corrector(Norm,H,dH);
550: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
551:
552: begin
553: Corrector(file,s,pivot,dist_sols,c,fail);
554: end Affine_Multiple_Severe_Normal_Reporting_Corrector;
555:
556: procedure Affine_Multiple_Severe_Conditioned_Silent_Corrector
557: ( s : in out Solu_Info_Array;
558: pivot : in out natural; dist_sols : in double_float;
559: c : in Corr_Pars; fail : out boolean ) is
560:
561: procedure Single_Corrector is
562: new Affine_Single_Severe_Conditioned_Silent_Corrector(Norm,H,dH);
563: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
564:
565: begin
566: Corrector(s,pivot,dist_sols,c,fail);
567: end Affine_Multiple_Severe_Conditioned_Silent_Corrector;
568:
569: procedure Affine_Multiple_Severe_Conditioned_Reporting_Corrector
570: ( file : in file_type; s : in out Solu_Info_Array;
571: pivot : in out natural; dist_sols : in double_float;
572: c : in Corr_Pars; fail : out boolean ) is
573:
574: procedure Single_Corrector is
575: new Affine_Single_Severe_Conditioned_Reporting_Corrector(Norm,H,dH);
576: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
577:
578: begin
579: Corrector(file,s,pivot,dist_sols,c,fail);
580: end Affine_Multiple_Severe_Conditioned_Reporting_Corrector;
581:
582: procedure Projective_Single_Loose_Normal_Silent_Corrector
583: ( s : in out Solu_Info; c : in Corr_Pars ) is
584:
585: info : integer;
586: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
587: j : Matrix(y'range,y'range);
588: ipvt : Standard_Natural_Vectors.Vector(y'range);
589: nit : natural := 0;
590: normv : double_float;
591:
592: begin
593: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
594: normv := Norm(s.sol.v);
595: if normv + 1.0 /= 1.0
596: then s.corr := s.cora / normv; s.resr := s.resa / normv;
597: end if;
598: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
599: j := dH(s.sol.v,s.sol.t);
600: for jj in y'range loop
601: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
602: end loop;
603: lufac(j,y'last,ipvt,info);
604: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
605: Min(y);
606: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
607: lusolve(j,y'last,ipvt,y); s.cora := Norm(y);
608: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
609: y := H(s.sol.v,s.sol.t);
610: y(y'last) := Create(0.0); s.resa := Norm(y);
611: normv := Norm(s.sol.v);
612: if normv + 1.0 /= 1.0
613: then s.corr := s.cora / normv; s.resr := s.resa / normv;
614: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
615: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
616: end loop;
617: end if;
618: nit := nit + 1;
619: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
620: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
621: end loop; -- STOP WHEN DESIRED PRECISION REACHED
622: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
623: exception
624: when numeric_error => return;
625: end Projective_Single_Loose_Normal_Silent_Corrector;
626:
627: procedure Projective_Single_Loose_Normal_Reporting_Corrector
628: ( file : in file_type;
629: s : in out Solu_Info; c : in Corr_Pars ) is
630:
631: info : integer;
632: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
633: j : Matrix(y'range,y'range);
634: ipvt : Standard_Natural_Vectors.Vector(y'range);
635: nit : natural := 0;
636: normv : double_float;
637:
638: begin
639: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
640: normv := Norm(s.sol.v);
641: if normv + 1.0 /= 1.0
642: then s.corr := s.cora / normv; s.resr := s.resa / normv;
643: end if;
644: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
645: j := dH(s.sol.v,s.sol.t);
646: for jj in y'range loop
647: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
648: end loop;
649: lufac(j,y'last,ipvt,info);
650: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
651: Min(y);
652: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
653: lusolve(j,y'last,ipvt,y); s.cora := Norm(y);
654: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
655: y := H(s.sol.v,s.sol.t);
656: y(y'last) := Create(0.0); s.resa := Norm(y);
657: normv := Norm(s.sol.v);
658: if normv + 1.0 /= 1.0
659: then s.corr := s.cora / normv; s.resr := s.resa / normv;
660: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
661: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
662: end loop;
663: end if;
664: cWrite(file,s.cora,s.corr,s.resa,s.resr); -- WRITE STATISTICS
665: nit := nit + 1;
666: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
667: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
668: end loop; -- STOP WHEN DESIRED PRECISION REACHED
669: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
670: exception
671: when numeric_error => return;
672: end Projective_Single_Loose_Normal_Reporting_Corrector;
673:
674: procedure Projective_Single_Loose_Conditioned_Silent_Corrector
675: ( s : in out Solu_Info; c : in Corr_Pars ) is
676:
677: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
678: n : constant natural := y'last;
679: j : Matrix(y'range,y'range);
680: qraux : Standard_Complex_Vectors.Vector(y'range)
681: := (y'range => Create(0.0));
682: dum : Standard_Complex_Vectors.Vector(y'range);
683: ipvt : Standard_Natural_Vectors.Vector(y'range);
684: info : integer;
685: nit : natural := 0;
686: normv : double_float;
687:
688: begin
689: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
690: normv := Norm(s.sol.v);
691: if normv + 1.0 /= 1.0
692: then s.corr := s.cora / normv; s.resr := s.resa / normv;
693: end if;
694: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
695: j := dH(s.sol.v,s.sol.t);
696: for jj in y'range loop
697: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
698: end loop;
699: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
700: Min(y); -- FOLLOWED BY LEAST SQUARES
701: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
702: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
703: s.cora := Norm(y);
704: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
705: y := H(s.sol.v,s.sol.t);
706: y(y'last) := Create(0.0); s.resa := Norm(y);
707: normv := Norm(s.sol.v);
708: if normv + 1.0 /= 1.0
709: then s.corr := s.cora / normv; s.resr := s.resa / normv;
710: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
711: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
712: end loop;
713: end if;
714: nit := nit + 1;
715: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
716: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
717: end loop; -- STOP WHEN DESIRED PRECISION REACHED
718: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
719: j := dH(s.sol.v,s.sol.t);
720: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
721: exception
722: when numeric_error => return;
723: end Projective_Single_Loose_Conditioned_Silent_Corrector;
724:
725: procedure Projective_Single_Loose_Conditioned_Reporting_Corrector
726: ( file : in file_type;
727: s : in out Solu_Info; c : in Corr_Pars ) is
728:
729: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
730: n : constant natural := y'last;
731: j : Matrix(y'range,y'range);
732: qraux : Standard_Complex_Vectors.Vector(y'range)
733: := (y'range => Create(0.0));
734: dum : Standard_Complex_Vectors.Vector(y'range);
735: ipvt : Standard_Natural_Vectors.Vector(y'range);
736: info : integer;
737: nit : natural := 0;
738: normv : double_float;
739:
740: begin
741: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
742: normv := Norm(s.sol.v);
743: if normv + 1.0 /= 1.0
744: then s.corr := s.cora / normv; s.resr := s.resa / normv;
745: end if;
746: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
747: j := dH(s.sol.v,s.sol.t);
748: for jj in y'range loop
749: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
750: end loop;
751: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
752: Min(y); -- FOLLOWED BY LEAST SQUARES
753: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
754: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
755: s.cora := Norm(y);
756: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
757: y := H(s.sol.v,s.sol.t);
758: y(y'last) := Create(0.0); s.resa := Norm(y);
759: normv := Norm(s.sol.v);
760: if normv + 1.0 /= 1.0
761: then s.corr := s.cora / normv; s.resr := s.resa / normv;
762: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
763: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
764: end loop;
765: end if;
766: cWrite(file,s.cora,s.corr,s.resa,s.resr); -- WRITE STATISTICS
767: cWrite(file,s.rcond,s.sol.m);
768: nit := nit + 1;
769: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
770: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
771: end loop; -- STOP WHEN DESIRED PRECISION REACHED
772: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
773: j := dH(s.sol.v,s.sol.t);
774: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
775: exception
776: when numeric_error => return;
777: end Projective_Single_Loose_Conditioned_Reporting_Corrector;
778:
779: procedure Projective_Single_Severe_Normal_Silent_Corrector
780: ( s : in out Solu_Info; c : in Corr_Pars ) is
781:
782: info : integer;
783: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
784: j : Matrix(y'range,y'range);
785: ipvt : Standard_Natural_Vectors.Vector(y'range);
786: nit : natural := 0;
787: normv,ncora,nresa,ncorr,nresr : double_float;
788: stop : boolean;
789:
790: begin
791: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
792: normv := Norm(s.sol.v);
793: if normv + 1.0 /= 1.0
794: then s.corr := s.cora / normv; s.resr := s.resa / normv;
795: end if;
796: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
797: j := dH(s.sol.v,s.sol.t);
798: for jj in y'range loop
799: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
800: end loop;
801: lufac(j,y'last,ipvt,info);
802: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
803: Min(y);
804: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
805: lusolve(j,y'last,ipvt,y); ncora := Norm(y);
806: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
807: y := H(s.sol.v,s.sol.t);
808: y(y'last) := Create(0.0); nresa := Norm(y);
809: normv := Norm(s.sol.v);
810: if normv + 1.0 /= 1.0
811: then ncorr := ncora / normv; nresr := nresa / normv;
812: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
813: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
814: end loop;
815: end if;
816: nit := nit + 1;
817: stop := (((ncora > s.cora) and then (nresa > s.resa))
818: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
819: -- STOP WHEN DIVERGENCE
820: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
821: exit when stop;
822: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
823: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
824: end loop; -- STOP WHEN DESIRED PRECISION REACHED
825: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
826: exception
827: when numeric_error => return;
828: end Projective_Single_Severe_Normal_Silent_Corrector;
829:
830: procedure Projective_Single_Severe_Normal_Reporting_Corrector
831: ( file : in file_type;
832: s : in out Solu_Info; c : in Corr_Pars ) is
833:
834: info : integer;
835: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
836: j : Matrix(y'range,y'range);
837: ipvt : Standard_Natural_Vectors.Vector(y'range);
838: nit : natural := 0;
839: normv,ncora,nresa,ncorr,nresr : double_float;
840: stop : boolean;
841:
842: begin
843: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
844: normv := Norm(s.sol.v);
845: if normv + 1.0 /= 1.0
846: then s.corr := s.cora / normv; s.resr := s.resa / normv;
847: end if;
848: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
849: j := dH(s.sol.v,s.sol.t);
850: for jj in y'range loop
851: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
852: end loop;
853: lufac(j,y'last,ipvt,info);
854: exit when info /= 0; -- STOP WHEN SINGULAR JACOBIAN
855: Min(y);
856: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
857: lusolve(j,y'last,ipvt,y); ncora := Norm(y);
858: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
859: y := H(s.sol.v,s.sol.t);
860: y(y'last) := Create(0.0); nresa := Norm(y);
861: normv := Norm(s.sol.v);
862: if normv + 1.0 /= 1.0
863: then ncorr := ncora / normv; nresr := nresa / normv;
864: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
865: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
866: end loop;
867: end if;
868: cWrite(file,ncora,ncorr,nresa,nresr); -- WRITE STATISTICS
869: nit := nit + 1;
870: stop := (((ncora > s.cora) and then (nresa > s.resa))
871: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
872: -- STOP WHEN DIVERGENCE
873: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
874: exit when stop;
875: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
876: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
877: end loop; -- STOP WHEN DESIRED PRECISION REACHED
878: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
879: exception
880: when numeric_error => return;
881: end Projective_Single_Severe_Normal_Reporting_Corrector;
882:
883: procedure Projective_Single_Severe_Conditioned_Silent_Corrector
884: ( s : in out Solu_Info; c : in Corr_Pars ) is
885:
886: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
887: n : constant natural := y'last;
888: j : Matrix(y'range,y'range);
889: qraux : Standard_Complex_Vectors.Vector(y'range)
890: := (y'range => Create(0.0));
891: dum : Standard_Complex_Vectors.Vector(y'range);
892: ipvt : Standard_Natural_Vectors.Vector(y'range);
893: info : integer;
894: nit : natural := 0;
895: normv,ncora,nresa,ncorr,nresr : double_float;
896: stop : boolean;
897:
898: begin
899: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
900: normv := Norm(s.sol.v);
901: if normv + 1.0 /= 1.0
902: then s.corr := s.cora / normv; s.resr := s.resa / normv;
903: end if;
904: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
905: j := dH(s.sol.v,s.sol.t);
906: for jj in y'range loop
907: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
908: end loop;
909: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
910: Min(y); -- FOLLOWED BY LEAST SQUARES
911: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
912: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
913: ncora := Norm(y);
914: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
915: y := H(s.sol.v,s.sol.t);
916: y(y'last) := Create(0.0); nresa := Norm(y);
917: normv := Norm(s.sol.v);
918: if normv + 1.0 /= 1.0
919: then ncorr := ncora / normv; nresr := nresa / normv;
920: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
921: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
922: end loop;
923: end if;
924: nit := nit + 1;
925: stop := (((ncora > s.cora) and then (nresa > s.resa))
926: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
927: -- STOP WHEN DIVERGENCE
928: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
929: exit when stop;
930: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
931: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
932: end loop; -- STOP WHEN DESIRED PRECISION REACHED
933: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
934: j := dH(s.sol.v,s.sol.t);
935: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
936: exception
937: when numeric_error => return;
938: end Projective_Single_Severe_Conditioned_Silent_Corrector;
939:
940: procedure Projective_Single_Severe_Conditioned_Reporting_Corrector
941: ( file : in file_type;
942: s : in out Solu_Info; c : in Corr_Pars ) is
943:
944: y : Vector(s.sol.v'range) := H(s.sol.v,s.sol.t);
945: n : constant natural := y'last;
946: j : Matrix(y'range,y'range);
947: qraux : Standard_Complex_Vectors.Vector(y'range)
948: := (y'range => Create(0.0));
949: dum : Standard_Complex_Vectors.Vector(y'range);
950: ipvt : Standard_Natural_Vectors.Vector(y'range);
951: info : integer;
952: nit : natural := 0;
953: normv,ncora,nresa,ncorr,nresr : double_float;
954: stop : boolean;
955:
956: begin
957: s.resa := Norm(y); s.cora := 1.0; -- INITIALIZATION
958: normv := Norm(s.sol.v);
959: if normv + 1.0 /= 1.0
960: then s.corr := s.cora / normv; s.resr := s.resa / normv;
961: end if;
962: while nit < c.maxit loop -- STOP WHEN MAXIMUM #ITERATIONS REACHED
963: j := dH(s.sol.v,s.sol.t);
964: for jj in y'range loop
965: j(j'last(1),jj) := s.sol.v(jj); -- CORRECT PERPENDICULAR
966: end loop;
967: QRD(j,qraux,ipvt,false); -- QR WITHOUT PIVOTING
968: Min(y); -- FOLLOWED BY LEAST SQUARES
969: y(y'last) := Create(0.0); -- IGNORE SCALING EQUATION
970: QRLS(j,n,n,n,qraux,y,dum,dum,y,dum,dum,100,info);
971: ncora := Norm(y);
972: Add(s.sol.v,y); s.length_path := s.length_path + s.cora;
973: y := H(s.sol.v,s.sol.t);
974: y(y'last) := Create(0.0); nresa := Norm(y);
975: normv := Norm(s.sol.v);
976: if normv + 1.0 /= 1.0
977: then ncorr := ncora / normv; nresr := nresa / normv;
978: for jj in s.sol.v'range loop -- SCALE THE SOLUTION
979: s.sol.v(jj) := s.sol.v(jj)/Create(normv);
980: end loop;
981: end if;
982: cWrite(file,ncora,ncorr,nresa,nresr); -- WRITE STATISTICS
983: cWrite(file,s.rcond,s.sol.m);
984: nit := nit + 1;
985: stop := (((ncora > s.cora) and then (nresa > s.resa))
986: or else ((ncorr > s.corr) and then (ncorr > s.corr)));
987: -- STOP WHEN DIVERGENCE
988: s.cora := ncora; s.resa := nresa; s.corr := ncorr; s.resr := nresr;
989: exit when stop;
990: exit when (((s.cora <= c.epsax) or else (s.resa <= c.epsaf))
991: and then ((s.corr <= c.epsrx) or else (s.resr <= c.epsrf)));
992: end loop; -- STOP WHEN DESIRED PRECISION REACHED
993: s.niter := s.niter + nit; s.nsyst := s.nsyst + nit;
994: j := dH(s.sol.v,s.sol.t);
995: lufco(j,y'last,ipvt,s.rcond); -- ESTIMATE CONDITION NUMBER
996: exception
997: when numeric_error => return;
998: end Projective_Single_Severe_Conditioned_Reporting_Corrector;
999:
1000: procedure Projective_Multiple_Loose_Normal_Silent_Corrector
1001: ( s : in out Solu_Info_Array;
1002: pivot : in out natural; dist_sols : in double_float;
1003: c : in Corr_Pars; fail : out boolean ) is
1004:
1005: procedure Single_Corrector is
1006: new Projective_Single_Loose_Normal_Silent_Corrector(Norm,H,dH);
1007: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
1008:
1009: begin
1010: Corrector(s,pivot,dist_sols,c,fail);
1011: end Projective_Multiple_Loose_Normal_Silent_Corrector;
1012:
1013: procedure Projective_Multiple_Loose_Normal_Reporting_Corrector
1014: ( file : in file_type; s : in out Solu_Info_Array;
1015: pivot : in out natural; dist_sols : in double_float;
1016: c : in Corr_Pars; fail : out boolean ) is
1017:
1018: procedure Single_Corrector is
1019: new Projective_Single_Loose_Normal_Reporting_Corrector(Norm,H,dH);
1020: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
1021:
1022: begin
1023: Corrector(file,s,pivot,dist_sols,c,fail);
1024: end Projective_Multiple_Loose_Normal_Reporting_Corrector;
1025:
1026: procedure Projective_Multiple_Loose_Conditioned_Silent_Corrector
1027: ( s : in out Solu_Info_Array;
1028: pivot : in out natural; dist_sols : in double_float;
1029: c : in Corr_Pars; fail : out boolean ) is
1030:
1031: procedure Single_Corrector is
1032: new Projective_Single_Loose_Conditioned_Silent_Corrector(Norm,H,dH);
1033: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
1034:
1035: begin
1036: Corrector(s,pivot,dist_sols,c,fail);
1037: end Projective_Multiple_Loose_Conditioned_Silent_Corrector;
1038:
1039: procedure Projective_Multiple_Loose_Conditioned_Reporting_Corrector
1040: ( file : in file_type; s : in out Solu_Info_Array;
1041: pivot : in out natural; dist_sols : in double_float;
1042: c : in Corr_Pars; fail : out boolean ) is
1043:
1044: procedure Single_Corrector is
1045: new Projective_Single_Loose_Normal_Reporting_Corrector(Norm,H,dH);
1046: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
1047:
1048: begin
1049: Corrector(file,s,pivot,dist_sols,c,fail);
1050: end Projective_Multiple_Loose_Conditioned_Reporting_Corrector;
1051:
1052: procedure Projective_Multiple_Severe_Normal_Silent_Corrector
1053: ( s : in out Solu_Info_Array;
1054: pivot : in out natural; dist_sols : in double_float;
1055: c : in Corr_Pars; fail : out boolean ) is
1056:
1057: procedure Single_Corrector is
1058: new Projective_Single_Severe_Normal_Silent_Corrector(Norm,H,dH);
1059: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
1060:
1061: begin
1062: Corrector(s,pivot,dist_sols,c,fail);
1063: end Projective_Multiple_Severe_Normal_Silent_Corrector;
1064:
1065: procedure Projective_Multiple_Severe_Normal_Reporting_Corrector
1066: ( file : in file_type; s : in out Solu_Info_Array;
1067: pivot : in out natural; dist_sols : in double_float;
1068: c : in Corr_Pars; fail : out boolean ) is
1069:
1070: procedure Single_Corrector is
1071: new Projective_Single_Loose_Normal_Reporting_Corrector(Norm,H,dH);
1072: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
1073:
1074: begin
1075: Corrector(file,s,pivot,dist_sols,c,fail);
1076: end Projective_Multiple_Severe_Normal_Reporting_Corrector;
1077:
1078: procedure Projective_Multiple_Severe_Conditioned_Silent_Corrector
1079: ( s : in out Solu_Info_Array;
1080: pivot : in out natural; dist_sols : in double_float;
1081: c : in Corr_Pars; fail : out boolean ) is
1082:
1083: procedure Single_Corrector is
1084: new Projective_Single_Severe_Conditioned_Silent_Corrector(Norm,H,dH);
1085: procedure Corrector is new Multiple_Silent_Corrector(Single_Corrector);
1086:
1087: begin
1088: Corrector(s,pivot,dist_sols,c,fail);
1089: end Projective_Multiple_Severe_Conditioned_Silent_Corrector;
1090:
1091: procedure Projective_Multiple_Severe_Conditioned_Reporting_Corrector
1092: ( file : in file_type; s : in out Solu_Info_Array;
1093: pivot : in out natural; dist_sols : in double_float;
1094: c : in Corr_Pars; fail : out boolean ) is
1095:
1096: procedure Single_Corrector is
1097: new Projective_Single_Severe_Conditioned_Reporting_Corrector(Norm,H,dH);
1098: procedure Corrector is new Multiple_Reporting_Corrector(Single_Corrector);
1099:
1100: begin
1101: Corrector(file,s,pivot,dist_sols,c,fail);
1102: end Projective_Multiple_Severe_Conditioned_Reporting_Corrector;
1103:
1104: end Correctors;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to correctors.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Continuation