![[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|
Return to homotopy.adb CVS log | Up to [local] / OpenXM_contrib / PHC / Ada / Homotopy |
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;