Annotation of OpenXM_contrib/PHC/Ada/Root_Counts/Implift/binomial_system_solvers.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Complex_Numbers; use Standard_Complex_Numbers;
! 2: with Standard_Complex_Numbers_Polar; use Standard_Complex_Numbers_Polar;
! 3: with Standard_Integer_Vectors;
! 4: with Standard_Complex_Norms_Equals; use Standard_Complex_Norms_Equals;
! 5: with Transforming_Solutions; use Transforming_Solutions;
! 6:
! 7: package body Binomial_System_Solvers is
! 8:
! 9: -- AN INTERMEDIATE ROUTINE : --
! 10:
! 11: procedure Factorize ( v : in out VecVec; n : in natural;
! 12: t : in out Transfo ) is
! 13:
! 14: -- DESCRIPTION :
! 15: -- This routines factorizes the binomial system defined by the
! 16: -- degrees in the Vector v;
! 17:
! 18: -- ON ENTRY :
! 19: -- v defines the degrees of binomial system
! 20: -- n the number of unknowns to be eliminated
! 21:
! 22: -- ON RETURN :
! 23: -- v the factorized binomial system
! 24: -- t the transformations used to factorize
! 25:
! 26: tt : Transfo;
! 27: pivot : integer;
! 28: tmp : Standard_Integer_Vectors.Link_to_Vector;
! 29:
! 30: begin
! 31: t := Id(v'last);
! 32: for i in 1..n loop
! 33: pivot := 0; -- search for pivot
! 34: for j in i..v'last loop
! 35: if v(j)(i) /= 0
! 36: then pivot := j;
! 37: exit;
! 38: end if;
! 39: end loop;
! 40: if pivot /= 0 -- interchange if necessary
! 41: then if pivot /= i
! 42: then tmp := v(i);
! 43: v(i) := v(pivot);
! 44: v(pivot) := tmp;
! 45: end if;
! 46: tt := Rotate(v(i),i);
! 47: Apply(tt,v(i));
! 48: for j in (i+1)..v'last loop
! 49: Apply(tt,v(j));
! 50: v(j).all(i) := 0;
! 51: end loop;
! 52: Mult1(t,tt);
! 53: Clear(tt);
! 54: end if;
! 55: end loop;
! 56: end Factorize;
! 57:
! 58: -- AUXILIARY ROUTINES :
! 59:
! 60: procedure Update ( sols1,sols2 : in out Solution_List ) is
! 61:
! 62: tmp : Solution_List := sols2;
! 63: ls : Link_to_Solution;
! 64:
! 65: begin
! 66: while not Is_Null(tmp) loop
! 67: ls := Head_Of(tmp);
! 68: declare
! 69: nls : Link_to_Solution := new Solution'(ls.all);
! 70: begin
! 71: Construct(nls,sols1);
! 72: end;
! 73: tmp := Tail_Of(tmp);
! 74: end loop;
! 75: Clear(sols2);
! 76: end Update;
! 77:
! 78: procedure Solve ( v : in VecVec; cv : in Vector;
! 79: i,n : in natural; sols : in out Solution_List;
! 80: acc : in out Vector ) is
! 81:
! 82: t : Transfo;
! 83: pivot,d : integer;
! 84: a : Complex_Number;
! 85: c : Vector(cv'range) := cv;
! 86: rv : Vector(cv'range);
! 87: tmp : Standard_Integer_Vectors.Link_to_Vector;
! 88: nv : VecVec(v'range);
! 89: temp : array(v'range) of integer;
! 90: newsols : Solution_List;
! 91:
! 92: begin
! 93: for j in i..v'last loop
! 94: nv(j) := new Standard_Integer_Vectors.Vector'(v(j).all);
! 95: end loop;
! 96: pivot := 0; -- search for pivot
! 97: for j in i..v'last loop
! 98: if nv(j)(i) /= 0
! 99: then pivot := j;
! 100: exit;
! 101: end if;
! 102: end loop;
! 103: if pivot /= 0 -- interchange if necessary
! 104: then if pivot /= i
! 105: then tmp := nv(i);
! 106: nv(i) := nv(pivot);
! 107: nv(pivot) := tmp;
! 108: a := c(i);
! 109: c(i) := c(pivot);
! 110: c(pivot) := a;
! 111: end if;
! 112: t := Rotate(nv(i),i);
! 113: Apply(t,nv(i));
! 114: for j in (i+1)..nv'last loop
! 115: Apply(t,nv(j));
! 116: temp(j) := nv(j)(i);
! 117: nv(j)(i) := 0;
! 118: end loop;
! 119: d := nv(i)(i); -- compute the solutions
! 120: if d < 0
! 121: then d := -d;
! 122: rv(i) := Create(1.0)/c(i);
! 123: else rv(i) := c(i);
! 124: end if;
! 125: for j in 1..d loop
! 126: a := Root(rv(i),d,j);
! 127: acc(i) := a;
! 128: for k in (i+1)..nv'last loop
! 129: rv(k) := c(k)*a**(-temp(k));
! 130: end loop;
! 131: if i < n
! 132: then Solve(nv,rv,i+1,n,newsols,acc);
! 133: else declare -- i = n
! 134: ls : Link_to_Solution;
! 135: s : Solution(n);
! 136: begin
! 137: s.t := Create(0.0);
! 138: s.m := 1;
! 139: s.v := acc;
! 140: ls := new Solution'(s);
! 141: Construct(ls,newsols);
! 142: end;
! 143: end if;
! 144: Transform(t,newsols);
! 145: Update(sols,newsols);
! 146: end loop;
! 147: Clear(t);
! 148: end if;
! 149: for j in i..nv'last loop
! 150: Standard_Integer_Vectors.Clear(nv(j));
! 151: end loop;
! 152: end Solve;
! 153:
! 154: -- THE SOLVER :
! 155:
! 156: procedure Solve ( v : in VecVec; cv : in Vector;
! 157: n : in natural; sols : in out Solution_List ) is
! 158:
! 159: wv : VecVec(v'range); -- workspace
! 160: acc : Vector(cv'range); -- accumulator
! 161:
! 162: begin
! 163: for i in v'range loop
! 164: wv(i) := new Standard_Integer_Vectors.Vector'(v(i).all);
! 165: end loop;
! 166: Solve(wv,cv,v'first,v'last,sols,acc);
! 167: Clear(wv);
! 168: end Solve;
! 169:
! 170: -- COMPUTING THE RESIDUALS :
! 171:
! 172: procedure Residuals ( v : in VecVec; cv : in Vector;
! 173: n : in natural; sols : in Solution_List;
! 174: res : out Vector ) is
! 175: nb : natural;
! 176: x : Vector(cv'range);
! 177: tmp : Solution_List;
! 178:
! 179: function Eval ( v : VecVec; cv,x : Vector ) return Vector is
! 180:
! 181: -- DESCRIPTION :
! 182: -- Computes the value of the vector x in the system defined by
! 183: -- v and cv.
! 184:
! 185: y : Standard_Complex_Vectors.Vector(x'range);
! 186:
! 187: begin
! 188: for i in x'range loop
! 189: y(i) := Create(1.0);
! 190: for j in v(i)'range loop
! 191: y(i) := y(i)*x(j)**v(i)(j);
! 192: end loop;
! 193: y(i) := y(i) - cv(i);
! 194: end loop;
! 195: return y;
! 196: end Eval;
! 197:
! 198: begin
! 199: tmp := sols;
! 200: nb := 1;
! 201: while not Is_Null(tmp) loop
! 202: x := Head_Of(tmp).v;
! 203: res(nb) := Create(Max_Norm(Eval(v,cv,x)));
! 204: tmp := Tail_Of(tmp);
! 205: nb := nb + 1;
! 206: end loop;
! 207: end Residuals;
! 208:
! 209: end Binomial_System_Solvers;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to binomial_system_solvers.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Root_Counts / Implift