Annotation of OpenXM_contrib/PHC/Ada/Homotopy/homotopy.adb, Revision 1.1
1.1 ! maekawa 1: with unchecked_deallocation;
! 2: with Standard_Floating_Numbers; use Standard_Floating_Numbers;
! 3: with Standard_Natural_Vectors;
! 4: with Standard_Complex_Vectors; use Standard_Complex_Vectors;
! 5: with Standard_Complex_Matrices; use Standard_Complex_Matrices;
! 6: with Standard_Complex_Polynomials; use Standard_Complex_Polynomials;
! 7: with Standard_Complex_Poly_Functions; use Standard_Complex_Poly_Functions;
! 8: with Standard_Complex_Poly_SysFun; use Standard_Complex_Poly_SysFun;
! 9: with Standard_Complex_Jaco_Matrices; use Standard_Complex_Jaco_Matrices;
! 10:
! 11: package body Homotopy is
! 12:
! 13: -- INTERNAL DATA -- to perform the numeric routines as fast as possible.
! 14:
! 15: type homtype is (nat,art); -- natural or artificial parameter homotopy
! 16:
! 17: type homdata ( ht : homtype; n,n1 : natural ) is record
! 18:
! 19: p : Poly_Sys(1..n); -- target system
! 20: pe : Eval_Poly_Sys(1..n); -- evaluable form of target
! 21: dh : Jaco_Mat(1..n,1..n1); -- Jacobian matrix of homotopy
! 22: dhe : Eval_Jaco_Mat(1..n,1..n1); -- evaluable form of dh
! 23:
! 24: case ht is
! 25: when nat =>
! 26: i : integer; -- which variable is parameter
! 27: when art =>
! 28: q,h : Poly_Sys(1..n); -- start system and homotopy
! 29: qe,he : Eval_Poly_Sys(1..n); -- evaluable form of q and h
! 30: dpe,dqe : Eval_Jaco_Mat(1..n,1..n); -- evaluable Jacobians
! 31: k : positive; -- relaxation power
! 32: gamma,beta : Vector(1..n); -- accessibility constants
! 33: linear : boolean; -- linear in t if true
! 34: end case;
! 35:
! 36: end record;
! 37:
! 38: type homtp is access homdata;
! 39: hom : homtp;
! 40:
! 41: -- GENERAL AUXILIARIES :
! 42:
! 43: function Mul ( v : Vector; p : Poly_Sys ) return Poly_Sys is
! 44:
! 45: -- DESCRIPTION :
! 46: -- Multiplies the ith polynomial with the ith entry in v.
! 47:
! 48: res : Poly_Sys(p'range);
! 49:
! 50: begin
! 51: for i in p'range loop
! 52: res(i) := v(i)*p(i);
! 53: end loop;
! 54: return res;
! 55: end Mul;
! 56:
! 57: function Mul ( m : Matrix; v : Vector ) return Matrix is
! 58:
! 59: -- DESCRIPTION :
! 60: -- Multiplies the ith row of m with the ith entry in v.
! 61:
! 62: res : Matrix(m'range(1),m'range(2));
! 63:
! 64: begin
! 65: for i in m'range(1) loop
! 66: for j in m'range(2) loop
! 67: res(i,j) := v(i)*m(i,j);
! 68: end loop;
! 69: end loop;
! 70: return res;
! 71: end Mul;
! 72:
! 73: -- AUXILIARIES TO THE CONSTRUCTORS :
! 74:
! 75: function Linear_Start_Factor
! 76: ( n,k : in natural; a : in Complex_Number ) return Poly is
! 77:
! 78: -- DESCRIPTION :
! 79: -- Returns a*(1-t)^k, with t as the (n+1)-th variable in the polynomial.
! 80:
! 81: res,tmp : Poly;
! 82: t : Term;
! 83:
! 84: begin
! 85: t.cf := a;
! 86: t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
! 87: res := Create(t); -- res = a
! 88: t.cf := Create(1.0);
! 89: tmp := Create(t); -- tmp = 1
! 90: t.dg(n+1) := 1;
! 91: Sub(tmp,t); -- tmp = 1-t
! 92: Standard_Natural_Vectors.Clear
! 93: (Standard_Natural_Vectors.Link_to_Vector(t.dg));
! 94: for i in 1..k loop
! 95: Mul(res,tmp);
! 96: end loop;
! 97: Clear(tmp);
! 98: return res;
! 99: end Linear_Start_Factor;
! 100:
! 101: function Nonlinear_Start_Factor
! 102: ( n,k : in natural; a : in Complex_Number ) return Poly is
! 103:
! 104: -- DESCRIPTION :
! 105: -- Returns (1 - (t - t*(1-t)*a))^k, with t as the (n+1)-th variable
! 106: -- in the polynomial.
! 107:
! 108: res,tmp : Poly;
! 109: t : Term;
! 110:
! 111: begin
! 112: t.cf := Create(1.0);
! 113: t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
! 114: res := Create(t); -- res = 1
! 115: tmp := Create(t); -- tmp = 1
! 116: t.dg(n+1) := 1;
! 117: Sub(tmp,t); -- tmp = 1-t
! 118: Mul(tmp,t); -- tmp = t*(1-t)
! 119: Mul(tmp,-a); -- tmp = -a*t*(1-t)
! 120: Add(tmp,t); -- tmp = t - a*t*(1-t)
! 121: Sub(res,tmp); -- res = 1 - t + a*t*(1-t)
! 122: if k > 1
! 123: then Copy(res,tmp);
! 124: for i in 1..(k-1) loop
! 125: Mul(res,tmp);
! 126: end loop;
! 127: end if;
! 128: Standard_Natural_Vectors.Clear
! 129: (Standard_Natural_Vectors.Link_to_Vector(t.dg));
! 130: Clear(tmp);
! 131: return res;
! 132: end Nonlinear_Start_Factor;
! 133:
! 134: function Linear_Target_Factor ( n,k : in natural ) return Poly is
! 135:
! 136: -- DESCRIPTION :
! 137: -- Returns t^k, with t as the (n+1)-th variable in the polynomial.
! 138:
! 139: res : Poly;
! 140: t : Term;
! 141:
! 142: begin
! 143: t.cf := Create(1.0);
! 144: t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
! 145: t.dg(n+1) := k;
! 146: res := Create(t);
! 147: Standard_Natural_Vectors.Clear
! 148: (Standard_Natural_Vectors.Link_to_Vector(t.dg));
! 149: return res;
! 150: end Linear_Target_Factor;
! 151:
! 152: function Linear_Target_Factor
! 153: ( n,k : in natural; a : in Complex_Number ) return Poly is
! 154:
! 155: -- DESCRIPTION :
! 156: -- Returns a*t^k, with t as the (n+1)-th variable in the polynomial.
! 157:
! 158: res : Poly;
! 159: t : Term;
! 160:
! 161: begin
! 162: t.cf := a;
! 163: t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
! 164: t.dg(n+1) := k;
! 165: res := Create(t);
! 166: Standard_Natural_Vectors.Clear
! 167: (Standard_Natural_Vectors.Link_to_Vector(t.dg));
! 168: return res;
! 169: end Linear_Target_Factor;
! 170:
! 171: function Nonlinear_Target_Factor
! 172: ( n,k : in natural; a : in Complex_Number ) return Poly is
! 173:
! 174: -- DESCRIPTION :
! 175: -- Returns (t - t*(1-t)*a)^k, with t as the (n+1)-th variable in
! 176: -- the polynomial.
! 177:
! 178: res,tmp : Poly;
! 179: t : Term;
! 180:
! 181: begin
! 182: t.cf := Create(1.0);
! 183: t.dg := new Standard_Natural_Vectors.Vector'(1..n+1 => 0);
! 184: tmp := Create(t); -- tmp = 1
! 185: t.dg(n+1) := 1;
! 186: res := Create(t); -- res = t
! 187: Sub(tmp,t); -- tmp = 1-t
! 188: Mul(tmp,t); -- tmp = t*(1-t)
! 189: Mul(tmp,a); -- tmp = a*t*(1-t)
! 190: Sub(res,tmp); -- res = t - a*t*(1-t)
! 191: if k > 1
! 192: then Copy(res,tmp);
! 193: for i in 1..(k-1) loop
! 194: Mul(res,tmp);
! 195: end loop;
! 196: end if;
! 197: Standard_Natural_Vectors.Clear
! 198: (Standard_Natural_Vectors.Link_to_Vector(t.dg));
! 199: Clear(tmp);
! 200: return res;
! 201: end Nonlinear_Target_Factor;
! 202:
! 203: function Plus_one_Unknown ( p : in Poly ) return Poly is
! 204:
! 205: -- DESCRIPTION :
! 206: -- The returning polynomial has place for one additional unknown.
! 207:
! 208: res : Poly;
! 209:
! 210: procedure Plus_Unknown_In_Term ( t : in out Term; c : out boolean ) is
! 211:
! 212: n : constant natural := t.dg'length;
! 213: temp : Standard_Natural_Vectors.Vector(1..(n+1));
! 214:
! 215: begin
! 216: temp(1..n) := t.dg.all;
! 217: temp(n+1) := 0;
! 218: Standard_Natural_Vectors.Clear
! 219: (Standard_Natural_Vectors.Link_to_Vector(t.dg));
! 220: t.dg := new Standard_Natural_Vectors.Vector'(temp);
! 221: c := true;
! 222: end Plus_Unknown_In_Term;
! 223: procedure Plus_Unknown_In_Terms is
! 224: new Changing_Iterator (process => Plus_Unknown_In_Term);
! 225:
! 226: begin
! 227: Copy(p,res);
! 228: Plus_Unknown_in_Terms(res);
! 229: return res;
! 230: end Plus_one_Unknown;
! 231:
! 232: function Linear_Homotopy
! 233: ( p,q : in Poly_Sys; k : in positive; a : in Complex_Number )
! 234: return Poly_Sys is
! 235:
! 236: h : Poly_Sys(p'range);
! 237: tempp,tempq,q_fac,p_fac : Poly;
! 238: n : natural := p'length;
! 239:
! 240: begin
! 241: q_fac := Linear_Start_Factor(n,k,a);
! 242: p_fac := Linear_Target_Factor(n,k);
! 243: for i in h'range loop
! 244: tempq := Plus_one_Unknown(q(i));
! 245: tempp := Plus_one_Unknown(p(i));
! 246: Mul(tempq,q_fac);
! 247: Mul(tempp,p_fac);
! 248: h(i) := tempq + tempp;
! 249: Clear(tempq); Clear(tempp);
! 250: end loop;
! 251: Clear(q_fac); Clear(p_fac);
! 252: return h;
! 253: end Linear_Homotopy;
! 254:
! 255: function Linear_Homotopy
! 256: ( p,q : in Poly_Sys; k : in positive; a : in Vector )
! 257: return Poly_Sys is
! 258:
! 259: h : Poly_Sys(p'range);
! 260: tempp,tempq,q_fac,p_fac : Poly;
! 261: n : natural := p'length;
! 262:
! 263: begin
! 264: for i in h'range loop
! 265: q_fac := Linear_Start_Factor(n,k,a(i));
! 266: p_fac := Linear_Target_Factor(n,k);
! 267: tempq := Plus_one_Unknown(q(i));
! 268: tempp := Plus_one_Unknown(p(i));
! 269: Mul(tempq,q_fac);
! 270: Mul(tempp,p_fac);
! 271: h(i) := tempq + tempp;
! 272: Clear(tempq); Clear(tempp);
! 273: Clear(q_fac); Clear(p_fac);
! 274: end loop;
! 275: return h;
! 276: end Linear_Homotopy;
! 277:
! 278: function Linear_Homotopy
! 279: ( p,q : in Poly_Sys; k : in positive; a,b : in Vector )
! 280: return Poly_Sys is
! 281:
! 282: h : Poly_Sys(p'range);
! 283: tempp,tempq,q_fac,p_fac : Poly;
! 284: n : natural := p'length;
! 285:
! 286: begin
! 287: for i in h'range loop
! 288: q_fac := Linear_Start_Factor(n,k,a(i));
! 289: p_fac := Linear_Target_Factor(n,k,b(i));
! 290: tempq := Plus_one_Unknown(q(i));
! 291: tempp := Plus_one_Unknown(p(i));
! 292: Mul(tempq,q_fac);
! 293: Mul(tempp,p_fac);
! 294: h(i) := tempq + tempp;
! 295: Clear(tempq); Clear(tempp);
! 296: Clear(q_fac); Clear(p_fac);
! 297: end loop;
! 298: return h;
! 299: end Linear_Homotopy;
! 300:
! 301: function Nonlinear_Homotopy
! 302: ( p,q : in Poly_Sys; k : in positive; a,b : in Vector )
! 303: return Poly_Sys is
! 304:
! 305: h : Poly_Sys(p'range);
! 306: tempp,tempq,q_fac,p_fac : Poly;
! 307: n : natural := p'length;
! 308:
! 309: begin
! 310: for i in h'range loop
! 311: q_fac := Nonlinear_Start_Factor(n,k,a(i));
! 312: p_fac := Nonlinear_Target_Factor(n,k,b(i));
! 313: tempq := Plus_one_Unknown(q(i));
! 314: tempp := Plus_one_Unknown(p(i));
! 315: Mul(tempq,q_fac);
! 316: Mul(tempp,p_fac);
! 317: h(i) := tempq + tempp;
! 318: Clear(tempq); Clear(tempp);
! 319: Clear(q_fac); Clear(p_fac);
! 320: end loop;
! 321: return h;
! 322: end Nonlinear_Homotopy;
! 323:
! 324: procedure Create ( p,q : in Poly_Sys; k : in positive;
! 325: a : in Complex_Number ) is
! 326:
! 327: n : constant natural := p'length;
! 328: dp, dq : Jaco_Mat(1..n,1..n);
! 329: ho : homdata(art,n,n+1);
! 330:
! 331: begin
! 332: Copy(p,ho.p); Copy(q,ho.q);
! 333: ho.h := Linear_Homotopy(p,q,k,a);
! 334: ho.pe := Create(ho.p);
! 335: ho.qe := Create(ho.q);
! 336: ho.he := Create(ho.h);
! 337: dp := Create(ho.p);
! 338: dq := Create(ho.q);
! 339: ho.dh := Create(ho.h);
! 340: ho.dpe := Create(dp);
! 341: ho.dqe := Create(dq);
! 342: ho.dhe := Create(ho.dh);
! 343: Clear(dp); Clear(dq);
! 344: ho.k := k;
! 345: for i in 1..n loop
! 346: ho.gamma(i) := a;
! 347: end loop;
! 348: for i in 1..n loop
! 349: ho.beta(i) := Create(1.0);
! 350: end loop;
! 351: ho.linear := true;
! 352: hom := new homdata'(ho);
! 353: end Create;
! 354:
! 355: procedure Create ( p,q : in Poly_Sys; k : in positive; a : in Vector ) is
! 356:
! 357: n : constant natural := p'length;
! 358: dp, dq : Jaco_Mat(1..n,1..n);
! 359: ho : homdata(art,n,n+1);
! 360:
! 361: begin
! 362: Copy(p,ho.p); Copy(q,ho.q);
! 363: ho.h := Linear_Homotopy(p,q,k,a);
! 364: ho.pe := Create(ho.p);
! 365: ho.qe := Create(ho.q);
! 366: ho.he := Create(ho.h);
! 367: dp := Create(ho.p);
! 368: dq := Create(ho.q);
! 369: ho.dh := Create(ho.h);
! 370: ho.dpe := Create(dp);
! 371: ho.dqe := Create(dq);
! 372: ho.dhe := Create(ho.dh);
! 373: Clear(dp); Clear(dq);
! 374: ho.k := k;
! 375: ho.gamma := a;
! 376: for i in 1..n loop
! 377: ho.beta(i) := Create(1.0);
! 378: end loop;
! 379: ho.linear := true;
! 380: hom := new homdata'(ho);
! 381: end Create;
! 382:
! 383: procedure Create ( p,q : in Poly_Sys; k : in positive; a,b : in Vector;
! 384: linear : in boolean ) is
! 385:
! 386: n : constant natural := p'length;
! 387: dp, dq : Jaco_Mat(1..n,1..n);
! 388: ho : homdata(art,n,n+1);
! 389:
! 390: begin
! 391: Copy(p,ho.p); Copy(q,ho.q);
! 392: ho.linear := linear;
! 393: if linear
! 394: then ho.h := Linear_Homotopy(p,q,k,a,b);
! 395: else ho.h := Nonlinear_Homotopy(p,q,k,a,b);
! 396: end if;
! 397: ho.pe := Create(ho.p);
! 398: ho.qe := Create(ho.q);
! 399: ho.he := Create(ho.h);
! 400: dp := Create(ho.p);
! 401: dq := Create(ho.q);
! 402: ho.dh := Create(ho.h);
! 403: ho.dpe := Create(dp);
! 404: ho.dqe := Create(dq);
! 405: ho.dhe := Create(ho.dh);
! 406: Clear(dp); Clear(dq);
! 407: ho.k := k;
! 408: ho.gamma := a;
! 409: ho.beta := b;
! 410: hom := new homdata'(ho);
! 411: end Create;
! 412:
! 413: procedure Create ( p : in Poly_Sys; k : in integer ) is
! 414:
! 415: n : constant natural := p'last-p'first+1;
! 416: ho : homdata(nat,n,n+1);
! 417:
! 418: begin
! 419: Copy(p,ho.p);
! 420: ho.pe := Create(ho.p);
! 421: ho.dh := Create(ho.p);
! 422: ho.dhe := Create(ho.dh);
! 423: ho.i := k;
! 424: hom := new homdata'(ho);
! 425: end Create;
! 426:
! 427: -- SELECTOR :
! 428:
! 429: function Homotopy_System return Poly_Sys is
! 430:
! 431: ho : homdata renames hom.all;
! 432:
! 433: begin
! 434: case ho.ht is
! 435: when nat => return ho.p;
! 436: when art => return ho.h;
! 437: end case;
! 438: end Homotopy_System;
! 439:
! 440: -- SYMBOLIC ROUTINES :
! 441:
! 442: function Eval ( t : Complex_Number ) return Poly_Sys is
! 443:
! 444: ho : homdata renames hom.all;
! 445:
! 446: begin
! 447: case ho.ht is
! 448: when nat => -- t = x(ho.i)
! 449: return Eval(ho.p,t,ho.i);
! 450: when art => -- compute : a * ((1 - t)^k) * q + (t^k) * p
! 451: declare
! 452: p_factor,q_factor : Vector(1..ho.n);
! 453: res,tmp : Poly_Sys(1..ho.n);
! 454: begin
! 455: if ho.linear
! 456: then
! 457: if AbsVal(t) = 0.0
! 458: then return Mul(ho.gamma,ho.q);
! 459: elsif abs(REAL_PART(t) - 1.0 ) + 1.0 = 1.0
! 460: and then abs(IMAG_PART(t)) + 1.0 = 1.0
! 461: then return Mul(ho.beta,ho.p);
! 462: else for i in 1..ho.n loop
! 463: q_factor(i) := ho.gamma(i);
! 464: p_factor(i) := ho.beta(i);
! 465: for i in 1..ho.k loop
! 466: q_factor(i) := (Create(1.0)-t) * q_factor(i);
! 467: p_factor(i) := t * p_factor(i);
! 468: end loop;
! 469: end loop;
! 470: res := Mul(p_factor,ho.p);
! 471: tmp := Mul(q_factor,ho.q);
! 472: Add(res,tmp);
! 473: Clear(tmp);
! 474: return res;
! 475: end if;
! 476: else return res; -- still to do !
! 477: end if;
! 478: end;
! 479: end case;
! 480: end Eval;
! 481:
! 482: function Diff ( t : Complex_Number ) return Poly_Sys is
! 483:
! 484: ho : homdata renames hom.all;
! 485:
! 486: begin
! 487: case ho.ht is
! 488: when nat => -- t = x(ho.i)
! 489: return Diff(ho.p,ho.i);
! 490: when art => -- compute - a*k*(1 - t)^(k-1)*q + b*k*t^(k-1)*p
! 491: declare
! 492: q_factor,p_factor : Vector(1..ho.n);
! 493: tmp : Poly_Sys(1..ho.n);
! 494: res : Poly_Sys(1..ho.n);
! 495: begin
! 496: if ho.linear
! 497: then
! 498: for i in 1..ho.n loop
! 499: q_factor(i) := (-ho.gamma(i)) * Create(double_float(ho.k));
! 500: p_factor(i) := ho.beta(i) * Create(double_float(ho.k));
! 501: end loop;
! 502: if AbsVal(t) = 0.0
! 503: then if ho.k = 1
! 504: then res := Mul(p_factor,ho.p);
! 505: tmp := Mul(q_factor,ho.q);
! 506: Add(res,tmp);
! 507: Clear(tmp);
! 508: return res;
! 509: else return Mul(q_factor,ho.q);
! 510: end if;
! 511: elsif abs( REAL_PART(t) - 1.0 ) + 1.0 = 1.0
! 512: and then abs(IMAG_PART(t)) + 1.0 = 1.0
! 513: then return Create(double_float(ho.k)) * ho.p;
! 514: else for i in 1..(ho.k-1) loop
! 515: q_factor := (Create(1.0)-t) * q_factor;
! 516: p_factor := t * p_factor;
! 517: end loop;
! 518: res := Mul(p_factor,ho.p);
! 519: tmp := Mul(q_factor,ho.q);
! 520: Add(res,tmp);
! 521: Clear(tmp);
! 522: return res;
! 523: end if;
! 524: else return res; -- still left to do !!!
! 525: end if;
! 526: end;
! 527: end case;
! 528: end Diff;
! 529:
! 530: -- NUMERIC ROUTINES :
! 531:
! 532: function Eval ( x : Vector; t : Complex_Number ) return Vector is
! 533:
! 534: ho : homdata renames hom.all;
! 535:
! 536: begin
! 537: case ho.ht is
! 538: when nat =>
! 539: declare
! 540: y : Vector(x'first..x'last+1);
! 541: begin
! 542: y(1..ho.i-1) := x(1..ho.i-1);
! 543: y(ho.i) := t;
! 544: y(ho.i+1..y'last) := x(ho.i..x'last);
! 545: return Eval(ho.pe,y);
! 546: end;
! 547: when art =>
! 548: if AbsVal(t) + 1.0 = 1.0
! 549: then if ho.linear
! 550: then return ho.gamma * Eval(ho.qe,x);
! 551: else return Eval(ho.qe,x);
! 552: end if;
! 553: elsif abs( REAL_PART(t) - 1.0 ) + 1.0 = 1.0
! 554: and then abs(IMAG_PART(t)) + 1.0 = 1.0
! 555: then if ho.linear
! 556: then return ho.beta * Eval(ho.pe,x);
! 557: else return Eval(ho.pe,x);
! 558: end if;
! 559: else declare
! 560: y : Vector(x'first..x'last+1);
! 561: begin
! 562: y(x'range) := x;
! 563: y(y'last) := t;
! 564: return Eval(ho.he,y);
! 565: end;
! 566: end if;
! 567: end case;
! 568: end Eval;
! 569:
! 570: function Diff ( x : Vector; t : Complex_Number ) return Vector is
! 571:
! 572: n : natural := x'length;
! 573:
! 574: begin
! 575: case hom.ht is
! 576: when nat => return Diff(x,t,hom.i);
! 577: when art => return Diff(x,t,n+1);
! 578: end case;
! 579: end Diff;
! 580:
! 581: function Diff ( x : Vector; t : Complex_Number ) return Matrix is
! 582:
! 583: ho : homdata renames hom.all;
! 584: n : natural renames ho.n;
! 585:
! 586: begin
! 587: case ho.ht is
! 588: when nat =>
! 589: declare
! 590: m : Matrix(1..n,1..n);
! 591: y : Vector(1..n+1);
! 592: begin
! 593: y(1..ho.i-1) := x(1..ho.i-1);
! 594: y(ho.i) := t;
! 595: y(ho.i+1..n+1) := x(ho.i..n);
! 596: for i in 1..n loop
! 597: for j in 1..n loop
! 598: m(i,j) := Eval(ho.dhe(i,j),y);
! 599: end loop;
! 600: end loop;
! 601: return m;
! 602: end;
! 603: when art =>
! 604: if AbsVal(t) + 1.0 = 1.0
! 605: then declare
! 606: m : Matrix(1..n,1..n) := Eval(ho.dqe,x);
! 607: begin
! 608: if ho.linear
! 609: then return Mul(m,ho.gamma);
! 610: else return m;
! 611: end if;
! 612: end;
! 613: elsif abs( REAL_PART(t) - 1.0 ) + 1.0 = 1.0
! 614: and then abs(IMAG_PART(t)) + 1.0 = 1.0
! 615: then declare
! 616: m : Matrix(1..n,1..n) := Eval(ho.dpe,x);
! 617: begin
! 618: if ho.linear
! 619: then return Mul(m,ho.beta);
! 620: else return m;
! 621: end if;
! 622: end;
! 623: else declare
! 624: m : Matrix(1..n,1..n);
! 625: y : Vector(1..n+1);
! 626: begin
! 627: y(1..n) := x;
! 628: y(n+1) := t;
! 629: for i in 1..n loop
! 630: for j in 1..n loop
! 631: m(i,j) := Eval(ho.dhe(i,j),y);
! 632: end loop;
! 633: end loop;
! 634: return m;
! 635: end;
! 636: end if;
! 637: end case;
! 638: end Diff;
! 639:
! 640: function Diff ( x : Vector; t : Complex_Number; k : natural ) return Vector is
! 641:
! 642: ho : homdata renames hom.all;
! 643: n : natural renames ho.n;
! 644: y : Vector(1..n+1);
! 645: res : Vector(1..n);
! 646:
! 647: begin
! 648: case ho.ht is
! 649: when nat => y(1..ho.k-1) := x(1..ho.i-1);
! 650: y(ho.i) := t;
! 651: y(ho.i+1..n+1) := x(ho.i..n);
! 652: when art => y(1..n) := x;
! 653: y(n+1) := t;
! 654: end case;
! 655: for i in 1..n loop
! 656: res(i) := Eval(ho.dhe(i,k),y);
! 657: end loop;
! 658: return res;
! 659: end Diff;
! 660:
! 661: function Diff ( x : Vector; t : Complex_Number; k : natural ) return Matrix is
! 662:
! 663: ho : homdata renames hom.all;
! 664: n : natural renames ho.n;
! 665: y : Vector(1..n+1);
! 666: res : Matrix(1..n,1..n);
! 667:
! 668: begin
! 669: case ho.ht is
! 670: when nat => y(1..ho.i-1) := x(1..ho.i-1);
! 671: y(ho.i) := t;
! 672: y(ho.i+1..n+1) := x(ho.i..n);
! 673: when art => y(1..n) := x;
! 674: y(n+1) := t;
! 675: end case;
! 676: for j in 1..(k-1) loop
! 677: for i in 1..n loop
! 678: res(i,j) := Eval(ho.dhe(i,j),y);
! 679: end loop;
! 680: end loop;
! 681: for j in (k+1)..(n+1) loop
! 682: for i in 1..n loop
! 683: res(i,j-1) := Eval(ho.dhe(i,j),y);
! 684: end loop;
! 685: end loop;
! 686: return res;
! 687: end Diff;
! 688:
! 689: -- DESTRUCTOR :
! 690:
! 691: procedure free is new unchecked_deallocation (homdata,homtp);
! 692:
! 693: procedure Clear is
! 694: begin
! 695: if hom /= null
! 696: then declare
! 697: ho : homdata renames hom.all;
! 698: begin
! 699: Clear(ho.p); Clear(ho.pe);
! 700: Clear(ho.dh); Clear(ho.dhe);
! 701: case ho.ht is
! 702: when nat => null;
! 703: when art =>
! 704: Clear(ho.q); Clear(ho.qe);
! 705: Clear(ho.h); Clear(ho.he);
! 706: Clear(ho.dpe); Clear(ho.dqe);
! 707: end case;
! 708: end;
! 709: free(hom);
! 710: end if;
! 711: end Clear;
! 712:
! 713: end Homotopy;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to homotopy.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Homotopy