Annotation of OpenXM_contrib/PHC/Ada/Continuation/directions_of_solution_paths.adb, Revision 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