Annotation of OpenXM_contrib/PHC/Ada/Continuation/directions_of_solution_paths.adb, Revision 1.1.1.1
1.1 maekawa 1: with integer_io; use integer_io;
2: with Standard_Floating_Numbers_io; use Standard_Floating_Numbers_io;
3: with Standard_Mathematical_Functions; use Standard_Mathematical_Functions;
4: with Standard_Floating_Vectors_io; use Standard_Floating_Vectors_io;
5: with Standard_Floating_Matrices; use Standard_Floating_Matrices;
6: with Standard_Floating_Matrices_io; use Standard_Floating_Matrices_io;
7: with Standard_Integer_Vectors;
8: with Standard_Integer_Vectors_io; use Standard_Integer_Vectors_io;
9: with vLpRs_Algorithm; use vLpRs_Algorithm;
10:
11: package body Directions_of_Solution_Paths is
12:
13: -- AUXILIARY :
14:
15: function Norm1 ( x : Standard_Floating_Vectors.Vector ) return double_float is
16:
17: res : double_float := 0.0;
18:
19: begin
20: for i in x'range loop
21: res := res + abs(x(i));
22: end loop;
23: return res;
24: end Norm1;
25:
26: procedure Shift_Up ( v : in out Standard_Floating_Vectors.Vector;
27: x : in double_float ) is
28:
29: -- DESCRIPTION :
30: -- Puts x at v(v'first) after a shift-up: v(i+1) := v(i).
31:
32: begin
33: for i in reverse v'first..(v'last-1) loop
34: v(i+1) := v(i);
35: end loop;
36: v(v'first) := x;
37: end Shift_Up;
38:
39: -- FIRST ORDER EXTRAPOLATION :
40:
41: procedure Affine_Update_Direction
42: ( t,prev_t,target : in Complex_Number;
43: x,prevx : in Standard_Complex_Vectors.Vector;
44: prevdls,prevstep : in out double_float;
45: prevdiff,v : in out Standard_Floating_Vectors.Vector ) is
46:
47: newdls : double_float;
48: newstep : double_float := AbsVal(t-prev_t);
49: newdiff : Standard_Floating_Vectors.Vector(prevdiff'range);
50: ratio,factor : double_float;
51:
52: begin
53: for i in v'range loop
54: newdiff(i) := LOG10(AbsVal(x(i))) - LOG10(AbsVal(prevx(i)));
55: end loop;
56: newdls := LOG10(AbsVal(target-prev_t)) - LOG10(AbsVal(target-t));
57: if prevdls /= 0.0
58: then ratio := prevstep/newstep;
59: factor := prevdls - ratio*newdls;
60: for i in v'range loop
61: v(i) := (prevdiff(i) - ratio*newdiff(i))/factor;
62: end loop;
63: end if;
64: prevdiff := newdiff;
65: prevstep := newstep;
66: prevdls := newdls;
67: end Affine_Update_Direction;
68:
69: procedure Projective_Update_Direction
70: ( t,prev_t,target : in Complex_Number;
71: x,prevx : in Standard_Complex_Vectors.Vector;
72: prevdls,prevstep : in out double_float;
73: prevdiff,v : in out Standard_Floating_Vectors.Vector ) is
74:
75: newdls : double_float;
76: newstep : double_float := AbsVal(t-prev_t);
77: newdiff : Standard_Floating_Vectors.Vector(prevdiff'range);
78: ratio,factor : double_float;
79:
80: begin
81: for i in v'range loop
82: newdiff(i) := ( LOG10(AbsVal(x(i))) - LOG10(AbsVal(x(x'last))) )
83: - ( LOG10(AbsVal(prevx(i))) - LOG10(AbsVal(prevx(prevx'last))) );
84: end loop;
85: newdls := LOG10(AbsVal(target-prev_t)) - LOG10(AbsVal(target-t));
86: if prevdls /= 0.0
87: then ratio := prevstep/newstep;
88: factor := prevdls - ratio*newdls;
89: for i in v'range loop
90: v(i) := (prevdiff(i) - ratio*newdiff(i))/factor;
91: end loop;
92: end if;
93: prevdiff := newdiff;
94: prevstep := newstep;
95: prevdls := newdls;
96: end Projective_Update_Direction;
97:
98: -- DATA MANAGEMENT FOR HIGHER-ORDER EXTRAPOLATION :
99:
100: procedure Extrapolation_Window
101: ( r,m : in natural; t,target : in Complex_Number;
102: x : in Standard_Complex_Vectors.Vector;
103: dt,s,logs : in out Standard_Floating_Vectors.Vector;
104: logx : in out VecVec ) is
105:
106: -- DESCRIPTION :
107: -- The data (dt,s,logs and logx) are stored as in a window: when
108: -- full at the incoming of a new element, all elements are shifted
109: -- and the oldest element drops off.
110:
111: use Standard_Floating_Vectors;
112:
113: begin
114: if (r = s'last) and (logx(r) /= null)
115: then for i in s'first+1..s'last loop -- shift the data
116: s(i-1) := s(i);
117: dt(i-1) := dt(i);
118: logs(i-1) := logs(i);
119: logx(i-1).all := logx(i).all;
120: end loop;
121: end if;
122: dt(r) := (AbsVal(target-t));
123: s(r) := (dt(r))**(1.0/double_float(m));
124: logs(r) := LOG10(s(r));
125: end Extrapolation_Window;
126:
127: procedure Refresh_Window ( r,m : in natural; dt : in Standard_Floating_Vectors.Vector;
128: s,logs : in out Standard_Floating_Vectors.Vector ) is
129:
130: -- DESCRIPTION :
131: -- Recomputes s and logs, after m has been changed.
132:
133: begin
134: for i in s'first..r loop
135: s(i) := (dt(i))**(1.0/double_float(m));
136: logs(i) := LOG10(s(i));
137: end loop;
138: end Refresh_Window;
139:
140: procedure Write_Update_Information
141: ( file : in file_type; r,m : in natural;
142: s,logs : in double_float;
143: logx : in Standard_Floating_Vectors.Vector ) is
144:
145: -- DESCRIPTION :
146: -- Writes r, s, log(s) and log|x(s)| on file.
147: -- The current version only writes a banner with r and m.
148:
149: begin
150: put(file,"extrapolation with r = "); put(file,r,1);
151: put(file," and m = "); put(file,m,1); --put_line(file," : ");
152: -- put(file,"m : "); put(file,s);
153: -- put(file," log(s) : "); put(file,logs);
154: -- put(file," log|x(s)| : ");
155: -- put(file,logx);
156: -- for i in logx'range loop
157: -- put(file,logx(i));
158: -- end loop;
159: new_line(file);
160: end Write_Update_Information;
161:
162: procedure Affine_Update_Extrapolation_Data
163: ( r,m : in natural; t,target : in Complex_Number;
164: x : in Standard_Complex_Vectors.Vector;
165: dt,s,logs : in out Standard_Floating_Vectors.Vector;
166: logx,wvl0,wvl1,wvltmp : in out VecVec ) is
167:
168: -- DESCRIPTION :
169: -- Updates the data needed to extrapolate in the affine case.
170:
171: use Standard_Floating_Vectors;
172:
173: begin
174: Extrapolation_Window(r,m,t,target,x,dt,s,logs,logx);
175: if logx(r) = null
176: then logx(r) := new Standard_Floating_Vectors.Vector(x'range);
177: end if;
178: if r > 0
179: then
180: if wvl0(r) = null
181: then wvl0(r) := new Standard_Floating_Vectors.Vector'(x'range => 0.0);
182: wvl1(r) := new Standard_Floating_Vectors.Vector'(x'range => 0.0);
183: wvltmp(r) := new Standard_Floating_Vectors.Vector'(x'range => 0.0);
184: end if;
185: end if;
186: for i in x'range loop
187: logx(r)(i) := LOG10(AbsVal(x(i)));
188: end loop;
189: end Affine_Update_Extrapolation_Data;
190:
191: procedure Projective_Update_Extrapolation_Data
192: ( r,m : in natural; t,target : in Complex_Number;
193: x : in Standard_Complex_Vectors.Vector;
194: dt,s,logs : in out Standard_Floating_Vectors.Vector;
195: logx : in out VecVec ) is
196:
197: -- DESCRIPTION :
198: -- Updates the data needed to extrapolate in the projective case,
199: -- under the assumption that the homogenization variable appears lastly.
200:
201: use Standard_Floating_Vectors;
202:
203: begin
204: Extrapolation_Window(r,m,t,target,x,dt,s,logs,logx);
205: if logx(r) = null
206: then logx(r) := new Standard_Floating_Vectors.Vector(x'first..x'last-1);
207: end if;
208: for i in x'first..x'last-1 loop
209: logx(r)(i) := LOG10(AbsVal(x(i))) - LOG10(AbsVal(x(x'last)));
210: end loop;
211: end Projective_Update_Extrapolation_Data;
212:
213: procedure Update_Errors
214: ( r : in natural;
215: errorv : in out Standard_Floating_Vectors.Vector;
216: error : out double_float; wvl0,wvl1,wvltmp : in out VecVec ) is
217:
218: -- DESCRIPTION :
219: -- Updates the error computation after the extrapolation.
220:
221: -- REQUIRED : r >= 1.
222:
223: begin
224: if r = 1
225: then for i in errorv'range loop
226: errorv(i) := abs(wvltmp(r)(i) - wvl1(r)(i));
227: end loop;
228: else for i in errorv'range loop
229: errorv(i) := abs(wvltmp(r-1)(i) - wvl1(r-1)(i));
230: end loop;
231: end if;
232: error := Norm1(errorv);
233: for i in 1..r loop
234: wvl0(i).all := wvl1(i).all;
235: wvl1(i).all := wvltmp(i).all;
236: end loop;
237: end Update_Errors;
238:
239: -- ESTIMATING THE CYCLE NUMBER :
240:
241: procedure Frequency_of_Estimate
242: ( newest,max : in natural; m,estm,cnt : in out natural;
243: eps : in double_float; newm : out boolean ) is
244:
245: -- DESCRIPTION :
246: -- This procedure manages the frequencies of the estimated values for m.
247: -- Only after the same estimate has been found a certain number of
248: -- times, the new estimate will be accepted.
249: -- The current version does not take the accuracy eps into account.
250:
251: -- ON ENTRY :
252: -- newest newly computed estimate for m;
253: -- max threshold on cnt before estm is returned;
254: -- m current value of m;
255: -- estm previous estimate;
256: -- cnt number of consecutive guesses that yielded estm;
257: -- eps accuracy of the current estimate.
258:
259: -- ON RETURN :
260: -- m new value of m;
261: -- estm new estimate;
262: -- cnt updated number of consecutive guesses that yielded estm;
263: -- newm true if m has changed.
264:
265: begin
266: if cnt = 0 -- initial estimate
267: then estm := newest;
268: cnt := 1;
269: elsif newest = estm -- update frequency for current estimate
270: then cnt := cnt + 1;
271: else cnt := 1; -- new estimate found
272: estm := newest;
273: end if;
274: if estm /= m -- look for modification of current cycle number
275: then if (cnt >= max) --and (eps <= 0.1)
276: then m := estm;
277: newm := true;
278: else newm := false;
279: end if;
280: else newm := false;
281: end if;
282: end Frequency_of_Estimate;
283:
284: procedure Extrapolate_on_Errors
285: ( file : in file_type;
286: r : in integer; h : in double_float;
287: err : in Standard_Floating_Vectors.Vector;
288: estm : out double_float ) is
289:
290: -- DESCRIPTION :
291: -- Performs an rth-order extrapolation on the errors.
292:
293: -- ON ENTRY :
294: -- file to write intermediate results on;
295: -- r order of the extrapolation method;
296: -- h ratio in geometric sequence;
297: -- err vector of range 0..r+1 with differences of estimates for
298: -- the outer normal, the most recent error is err(0).
299:
300: -- ON RETURN :
301: -- extm estimated value for m.
302:
303: em,hm,exterr : Standard_Floating_Vectors.Vector(1..r+1);
304: dlog : constant double_float := log10(h);
305: f : double_float;
306:
307: begin
308: for j in exterr'range loop
309: exterr(j) := log10(err(j-1)) - log10(err(j));
310: end loop;
311: em(1) := dlog/exterr(1); -- 0th order estimate
312: -- if (m(1) < 0.0001) or (m(1) > 1000.0) -- avoid divergence
313: -- then m(r+1) := m(1);
314: -- else
315: hm(1) := h**(1.0/em(1));
316: for k in 1..r loop
317: f := hm(k) - 1.0;
318: for j in 1..r-k+1 loop
319: exterr(j) := exterr(j+1) + (exterr(j+1) - exterr(j))/f;
320: end loop;
321: em(k+1) := dlog/exterr(1);
322: -- exit when ((m(k+1) < 0.0001) or (m(k+1) > 1000.0));
323: hm(k+1) := h**(double_float(k+1)/em(k+1));
324: end loop;
325: -- end if;
326: estm := em(r+1);
327: put(file,"em(0.."); put(file,r,1); put(file,") : ");
328: for i in em'range loop
329: put(file,em(i),3,3,3);
330: end loop;
331: new_line(file);
332: exception
333: when others => null;
334: end Extrapolate_on_Errors;
335:
336: procedure Estimate0
337: ( r,max : in natural; m,estm,cnt : in out natural;
338: dt : in Standard_Floating_Vectors.Vector;
339: err,newerr : in double_float; rat,eps : out double_float;
340: newm : out boolean ) is
341:
342: -- DESCRIPTION :
343: -- Returns an 0th-order estimate of the cycle number m.
344:
345: -- ON ENTRY :
346: -- r extrapolation order;
347: -- max threshold on cnt before estm is returned;
348: -- m current value of m;
349: -- estm previous estimate;
350: -- cnt number of consecutive guesses that yielded estm;
351: -- dt consecutive distances to the target;
352: -- err previous error;
353: -- newerr current error.
354:
355: -- ON RETURN :
356: -- m new value of m;
357: -- estm new estimate;
358: -- cnt updated number of consecutive guesses that yielded estm;
359: -- rat ratio used to estimate m;
360: -- eps accuracy of the rounding value for m;
361: -- newm true if m has changed.
362:
363: dferr : constant double_float := log10(newerr) - log10(err);
364: h : constant double_float := dt(r)/dt(r-1);
365: dfsr1 : constant double_float := log10(h);
366: ratio : constant double_float := abs(dfsr1/dferr);
367: res : natural := integer(ratio);
368: accuracy : constant double_float := abs(double_float(res) - ratio);
369:
370: begin
371: if res = 0
372: then res := 1;
373: end if;
374: Frequency_of_Estimate(res,max,m,estm,cnt,accuracy,newm);
375: rat := ratio; eps := accuracy;
376: end Estimate0;
377:
378: procedure Estimate
379: ( file : in file_type; r : in integer;
380: max : in natural; m,estm,cnt : in out natural;
381: h : in double_float;
382: diferr : in Standard_Floating_Vectors.Vector;
383: rat,eps : out double_float; newm : out boolean ) is
384:
385: -- DESCRIPTION :
386: -- Estimates m by extrapolation on consecutvie differences of errors,
387: -- stored in the parameter diferr. For the specfication of the other
388: -- parameters, see the procedure Estimate0.
389:
390: res : integer;
391: fltestm,accuracy : double_float;
392:
393: begin
394: -- if r < dt'last
395: -- then Extrapolate_on_Errors(file,r-1,h,diferr(0..r),fltestm);
396: -- else
397: Extrapolate_on_Errors(file,r,h,diferr(0..r+1),fltestm);
398: -- end if;
399: res := integer(fltestm);
400: if res <= 0
401: then res := 1;
402: end if;
403: accuracy := abs(double_float(res) - fltestm);
404: Frequency_of_Estimate(res,max,m,estm,cnt,accuracy,newm);
405: rat := fltestm; eps := accuracy;
406: end Estimate;
407:
408: -- APPLYING THE vLpRs-Algorithm :
409:
410: function vLpRs_Extrapolate
411: ( r : natural; s,logs : Standard_Floating_Vectors.Vector;
412: logx : VecVec )
413: return Standard_Floating_Vectors.Vector is
414:
415: -- DESCRIPTION :
416: -- Higher-order extrapolation is based on vLpRs-Algorithm.
417:
418: -- REQUIRED : r >= 1.
419:
420: res : Standard_Floating_Vectors.Vector(logx(r).all'range);
421: logx1 : Standard_Floating_Vectors.Vector(0..r);
422: rt1,rt2 : Matrix(1..r-1,1..r-1);
423: srp,dsp : Standard_Floating_Vectors.Vector(1..r-1) := (0..r-1 => 0.0);
424: p : Standard_Floating_Vectors.Vector(0..r-1) := (0..r-1 => 0.0);
425: l,v : Standard_Floating_Vectors.Vector(0..r) := (0..r => 1.0);
426:
427: begin
428: rt1(1,1) := 0.0; rt2(1,1) := 0.0;
429: for i in res'range loop
430: for j in logx1'range loop
431: logx1(j) := logx(j)(i);
432: end loop;
433: vLpRs_full(r,s(0..r),logs(0..r),logx1(0..r),srp,dsp,p,l,v,rt1,rt2);
434: res(i) := v(r)/l(r);
435: end loop;
436: return res;
437: end vLpRs_Extrapolate;
438:
439: function vLpRs_Extrapolate
440: ( file : file_type; r : natural;
441: s,logs : Standard_Floating_Vectors.Vector; logx : VecVec )
442: return Standard_Floating_Vectors.Vector is
443:
444: -- DESCRIPTION :
445: -- Higher-order extrapolation based on vLpRs-Algorithm.
446:
447: -- REQUIRED : r >= 1.
448:
449: res : Standard_Floating_Vectors.Vector(logx(r).all'range);
450: logx1 : Standard_Floating_Vectors.Vector(0..r);
451: rt1,rt2 : Matrix(1..r-1,1..r-1);
452: srp,dsp : Standard_Floating_Vectors.Vector(1..r-1) := (0..r-1 => 0.0);
453: p : Standard_Floating_Vectors.Vector(0..r-1) := (0..r-1 => 0.0);
454: l,v : Standard_Floating_Vectors.Vector(0..r) := (0..r => 0.0);
455:
456: begin
457: for i in rt1'range(1) loop
458: for j in rt1'range(2) loop
459: rt1(i,j) := 0.0; rt2(i,j) := 0.0;
460: end loop;
461: end loop;
462: for i in res'range loop
463: for j in logx1'range loop
464: logx1(j) := logx(j)(i);
465: end loop;
466: vLpRs_pipe(file,r,s(0..r),logs(0..r),logx1(0..r),srp,dsp,p,l,v,rt1,rt2);
467: res(i) := v(r)/l(r);
468: end loop;
469: return res;
470: end vLpRs_Extrapolate;
471:
472: procedure vLpRs_Extrapolate
473: ( file : in file_type; r : in natural;
474: s,logs : in Standard_Floating_Vectors.Vector;
475: logx,wvl0 : in VecVec; wvl1 : in out VecVec;
476: w,wv,wl : out Standard_Floating_Vectors.Vector ) is
477:
478: -- DESCRIPTION :
479: -- Higher-order extrapolation based on vLpRs-Algorithm,
480: -- writes an error table on file.
481:
482: -- REQUIRED : r >= 1.
483:
484: n : constant natural := logx(r)'length;
485: logx1 : Standard_Floating_Vectors.Vector(0..r);
486: rt1,rt2 : Matrix(1..r-1,1..r-1);
487: srp,dsp : Standard_Floating_Vectors.Vector(1..r-1) := (1..r-1 => 0.0);
488: p : Standard_Floating_Vectors.Vector(0..r-1) := (0..r-1 => 0.0);
489: l,v : Standard_Floating_Vectors.Vector(0..r) := (0..r => 0.0);
490: error : Standard_Floating_Vectors.Vector(1..r) := (1..r => 0.0);
491:
492: begin
493: for i in rt1'range(1) loop
494: for j in rt1'range(2) loop
495: rt1(i,j) := 0.0; rt2(i,j) := 0.0;
496: end loop;
497: end loop;
498: for i in logx(r)'range loop
499: for j in logx1'range loop
500: logx1(j) := logx(j)(i);
501: end loop;
502: vLpRs_pipe(file,r,s(0..r),logs(0..r),logx1(0..r),srp,dsp,p,l,v,rt1,rt2);
503: w(i) := v(r)/l(r);
504: put(file,s(r),2,3,3);
505: for j in 1..r loop
506: wvl1(j)(i) := v(j)/l(j);
507: error(j) := abs(wvl1(j)(i)-wvl0(j)(i));
508: put(file,error(j),2,3,3);
509: end loop;
510: new_line(file);
511: end loop;
512: wv := v; wl := l;
513: end vLpRs_Extrapolate;
514:
515: -- HIGHER-ORDER EXTRAPOLATION (target routines) :
516:
517: procedure Affine_Update_Direction
518: ( r,m,estm,cntm : in out natural; thresm : in natural;
519: er : in out integer; t,target : in Complex_Number;
520: x : in Standard_Complex_Vectors.Vector;
521: dt,s,logs : in out Standard_Floating_Vectors.Vector;
522: logx,wvl0,wvl1,wvltmp : in out VecVec;
523: v,diferr : in out Standard_Floating_Vectors.Vector;
524: error : in out double_float ) is
525:
526: use Standard_Floating_Vectors;
527: res : Standard_Floating_Vectors.Vector(v'range);
528: eps,newerr : double_float;
529: newm : boolean := false;
530:
531: begin
532: Affine_Update_Extrapolation_Data
533: (r,m,t,target,x,dt,s,logs,logx,wvl0,wvl1,wvltmp);
534: if r >= 1
535: then res := vLpRs_Extrapolate(r,s,logs,logx);
536: newerr := Norm1(res-v);
537: end if;
538: if r < s'last
539: then r := r+1;
540: end if;
541: if r >= 3 and (newerr < error)
542: then --Estimate(r,r,thresm,m,estm,cntm,dt,s,logs,error,newerr,eps,newm);
543: if newm
544: then res := vLpRs_Extrapolate(r,s,logs,logx);
545: newerr := Norm1(res-v);
546: end if;
547: end if;
548: if r >= 1
549: then v := res;
550: error := newerr;
551: end if;
552: end Affine_Update_Direction;
553:
554: procedure Affine_Update_Direction
555: ( file : in file_type;
556: r,m,estm,cntm : in out natural; thresm : in natural;
557: er : in out integer; t,target : in Complex_Number;
558: x : in Standard_Complex_Vectors.Vector;
559: dt,s,logs : in out Standard_Floating_Vectors.Vector;
560: logx,wvl0,wvl1,wvltmp : in out VecVec;
561: v,diferr : in out Standard_Floating_Vectors.Vector;
562: error : in out double_float ) is
563:
564: use Standard_Floating_Vectors;
565: res,errorv,newerrv : Standard_Floating_Vectors.Vector(v'range);
566: wv,wl : Standard_Floating_Vectors.Vector(0..r);
567: ind : natural := 1; -- errors based on first-order extrapolation
568: rat,eps,newerr : double_float;
569: mv,cntmv,estmv : Standard_Integer_Vectors.Vector(v'range);
570: newm : boolean := false;
571:
572: begin
573: Affine_Update_Extrapolation_Data
574: (r,m,t,target,x,dt,s,logs,logx,wvl0,wvl1,wvltmp);
575: Write_Update_Information(file,r,m,s(r),logs(r),logx(r).all);
576: if r >= 1
577: then vLpRs_Extrapolate(file,r,s,logs,logx,wvl1,wvltmp,res,wv,wl);
578: if r = 1 then diferr(0) := 1.0; end if;
579: for i in errorv'range loop
580: errorv(i) := abs(wvltmp(ind)(i) - wvl1(ind)(i));
581: end loop;
582: Shift_Up(diferr,Norm1(errorv));
583: if er < diferr'last-1 then er := er+1; end if;
584: end if;
585: if er >= 1 and (diferr(0) < diferr(1))
586: then-- Estimate0(r,thresm,m,estm,cntm,diferr(1),diferr(0),rat,eps,newm);
587: Estimate(file,er,thresm,m,estm,cntm,dt(r)/dt(r-1),
588: diferr,rat,eps,newm);
589: put(file,"Ratio for m : "); put(file,rat,3,3,3);
590: put(file," and accuracy : "); put(file,eps,3,3,3); new_line(file);
591: if newm
592: then Refresh_Window(r,m,dt,s,logs);
593: put_line(file,"Extrapolation after adjusting the m-value :");
594: vLpRs_Extrapolate(file,r,s,logs,logx,wvl1,wvltmp,res,wv,wl);
595: for i in errorv'range loop
596: errorv(i) := abs(wvltmp(ind)(i) - wvl1(ind)(i));
597: end loop;
598: Shift_Up(diferr,Norm1(errorv)); er := -2;
599: put(file,"old direction : "); put(file,v); new_line(file);
600: put(file,"new direction : "); put(file,res); new_line(file);
601: end if;
602: end if;
603: if r >= 1
604: then v := res;
605: Update_Errors(r,errorv,error,wvl0,wvl1,wvltmp);
606: end if;
607: if r < s'last then r := r+1; end if;
608: end Affine_Update_Direction;
609:
610: procedure Projective_Update_Direction
611: ( r,m,estm,cntm : in out natural; thresm : in natural;
612: er : in out integer; t,target : in Complex_Number;
613: x : in Standard_Complex_Vectors.Vector;
614: dt,s,logs : in out Standard_Floating_Vectors.Vector;
615: logx : in out VecVec;
616: prevv,v : in out Standard_Floating_Vectors.Vector;
617: error : in out double_float ) is
618:
619: use Standard_Floating_Vectors;
620: res : Standard_Floating_Vectors.Vector(v'range);
621: eps,newerr : double_float;
622: newm : boolean := false;
623:
624: begin
625: Projective_Update_Extrapolation_Data(r,m,t,target,x,dt,s,logs,logx);
626: if r >= 1
627: then res := vLpRs_Extrapolate(r,s,logs,logx);
628: error := Norm1(res-prevv); newerr := Norm1(res-v);
629: prevv := v;
630: v := res;
631: end if;
632: -- if r >= 2 and (newerr < error)
633: -- then Estimate(r,r,thresm,m,estm,cntm,dt,s,logs,error,newerr,eps,newm);
634: -- end if;
635: if r < s'last
636: then r := r+1;
637: end if;
638: error := newerr;
639: end Projective_Update_Direction;
640:
641: procedure Projective_Update_Direction
642: ( file : in file_type;
643: r,m,estm,cntm : in out natural; thresm : in natural;
644: er : in out integer; t,target : in Complex_Number;
645: x : in Standard_Complex_Vectors.Vector;
646: dt,s,logs : in out Standard_Floating_Vectors.Vector;
647: logx : in out VecVec;
648: prevv,v : in out Standard_Floating_Vectors.Vector;
649: error : in out double_float ) is
650:
651: use Standard_Floating_Vectors;
652: res : Standard_Floating_Vectors.Vector(v'range);
653: eps,newerr : double_float;
654: newm : boolean := false;
655:
656: begin
657: Projective_Update_Extrapolation_Data(r,m,t,target,x,dt,s,logs,logx);
658: Write_Update_Information(file,r,m,s(r),logs(r),logx(r).all);
659: if r >= 1
660: then res := vLpRs_Extrapolate(file,r,s,logs,logx);
661: error := Norm1(res-prevv);
662: newerr := Norm1(res-v);
663: prevv := v;
664: v := res;
665: end if;
666: if r < s'last
667: then r := r+1;
668: end if;
669: -- if r >= 2 and (newerr < error)
670: -- then Estimate(r,r,thresm,m,estm,cntm,dt,s,logs,error,newerr,eps,newm);
671: -- end if;
672: error := newerr;
673: end Projective_Update_Direction;
674:
675: end Directions_of_Solution_Paths;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to directions_of_solution_paths.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