Annotation of OpenXM_contrib/PHC/Ada/Math_Lib/Supports/face_enumerators_utilities.adb, Revision 1.1
1.1 ! maekawa 1: with Standard_Common_Divisors; use Standard_Common_Divisors;
! 2: with Standard_Integer_Matrices; use Standard_Integer_Matrices;
! 3: with Standard_Integer_Linear_Solvers; use Standard_Integer_Linear_Solvers;
! 4:
! 5: package body Face_Enumerators_Utilities is
! 6:
! 7: function Is_Zero ( v : Vector ) return boolean is
! 8: begin
! 9: for i in v'range loop
! 10: if v(i) /= 0
! 11: then return false;
! 12: end if;
! 13: end loop;
! 14: return true;
! 15: end Is_Zero;
! 16:
! 17: function gcd ( v : Vector ) return integer is
! 18:
! 19: g : integer := v(v'first);
! 20:
! 21: begin
! 22: if g < 0
! 23: then g := -g;
! 24: end if;
! 25: if g /= 1
! 26: then for i in v'first+1..v'last loop
! 27: if v(i) /= 0
! 28: then if v(i) < 0
! 29: then g := gcd(g,-v(i));
! 30: else g := gcd(g,v(i));
! 31: end if;
! 32: end if;
! 33: exit when (g = 1);
! 34: end loop;
! 35: end if;
! 36: return g;
! 37: end gcd;
! 38:
! 39: procedure Scale ( v : in out Vector ) is
! 40:
! 41: g : integer := gcd(v);
! 42:
! 43: begin
! 44: if (g /= 0) and then (g /= 1)
! 45: then for i in v'range loop
! 46: v(i) := v(i)/g;
! 47: end loop;
! 48: end if;
! 49: end Scale;
! 50:
! 51: function Is_In ( x : integer; v : Vector ) return boolean is
! 52: begin
! 53: for k in v'range loop
! 54: if v(k) = x
! 55: then return true;
! 56: end if;
! 57: end loop;
! 58: return false;
! 59: end Is_In;
! 60:
! 61: function In_Edge ( x,a,b : Vector ) return boolean is
! 62:
! 63: ba,xa,piv : integer;
! 64:
! 65: begin
! 66: for i in b'range loop -- look for first i: b(i) - a(i) /= 0
! 67: ba := b(i) - a(i);
! 68: if ba /= 0
! 69: then piv := i; exit;
! 70: end if;
! 71: end loop;
! 72: if ba = 0
! 73: then return Equal(x,a); -- in this case b = a, so test x = a
! 74: else if ba < 0
! 75: then ba := -ba;
! 76: xa := x(piv) - b(piv);
! 77: else xa := x(piv) - a(piv);
! 78: end if;
! 79: if xa*ba >= 0 and then xa <= ba
! 80: then return true;
! 81: else return false;
! 82: end if;
! 83: end if;
! 84: end In_Edge;
! 85:
! 86: function Last_Rows_Zero ( mat : matrix; frst : integer ) return boolean is
! 87:
! 88: -- DESCRIPTION :
! 89: -- Returns true when the last rows of mat, starting at the index frst,
! 90: -- are completely zero.
! 91:
! 92: function Zero_Row ( mat : matrix; i : integer ) return boolean is
! 93: begin
! 94: for j in mat'range(2) loop
! 95: if mat(i,j) /= 0
! 96: then return false;
! 97: end if;
! 98: end loop;
! 99: return true;
! 100: end Zero_Row;
! 101:
! 102: begin
! 103: for i in frst..mat'last(1) loop
! 104: if not Zero_Row(mat,i)
! 105: then return false;
! 106: end if;
! 107: end loop;
! 108: return true;
! 109: end Last_Rows_Zero;
! 110:
! 111: function In_Face ( k : in natural; face,x : Vector; pts : VecVec )
! 112: return boolean is
! 113:
! 114: mat : Matrix(x'range,1..k+1);
! 115: sol : Vector(1..k+1) := (1..k+1 => 0);
! 116: sum : integer := 0;
! 117:
! 118: begin
! 119: if k = 0
! 120: then return (pts(face(face'first)).all = x);
! 121: elsif k = 1
! 122: then return In_Edge(x,pts(face(face'first)).all,
! 123: pts(face(face'first+1)).all);
! 124: end if;
! 125: -- put_line("Calling In_Face for ");
! 126: -- put("x : "); put(x); new_line;
! 127: -- for i in face'range loop
! 128: -- put("pts("); put(face(i),1); put(") : "); put(pts(face(i))); new_line;
! 129: -- end loop;
! 130: -- compute the decomposition of x w.r.t. the face :
! 131: for j in 1..k loop
! 132: for i in mat'range(1) loop
! 133: mat(i,j) := pts(face(j))(i) - pts(face(0))(i);
! 134: end loop;
! 135: end loop;
! 136: for i in mat'range(1) loop
! 137: mat(i,k+1) := x(i) - pts(face(0))(i);
! 138: end loop;
! 139: -- put_line("The matrix which defines the decomposition :");
! 140: Upper_Triangulate(mat);
! 141: -- put(mat);
! 142: if not Last_Rows_Zero(mat,k+2)
! 143: then return false;
! 144: else declare
! 145: matk1 : matrix(1..k+1,1..k+1);
! 146: begin
! 147: for i in matk1'range(1) loop
! 148: for j in matk1'range(2) loop
! 149: matk1(i,j) := mat(i,j);
! 150: end loop;
! 151: end loop;
! 152: Scale(matk1); Solve0(matk1,sol);
! 153: -- put("The computed solution : "); put(sol); new_line;
! 154: for i in 1..k loop
! 155: if sol(i)*sol(k+1) > 0
! 156: then --put_line("x lies not in the face");
! 157: return false;
! 158: else sum := sum + sol(i);
! 159: end if;
! 160: end loop;
! 161: -- put("The sum : "); put(sum,1); new_line;
! 162: if sum = 0
! 163: then --put_line("x lies not in the face");
! 164: return false;
! 165: end if;
! 166: if sol(k+1) < 0
! 167: then sol(k+1) := -sol(k+1);
! 168: else sum := -sum;
! 169: end if;
! 170: if sum > sol(k+1)
! 171: then --put_line("x lies not in the face");
! 172: return false;
! 173: else --put_line("x lies in the face");
! 174: return true;
! 175: end if;
! 176: end;
! 177: end if;
! 178: end In_Face;
! 179:
! 180: end Face_Enumerators_Utilities;
FreeBSD-CVSweb <freebsd-cvsweb@FreeBSD.org>
![[BACK]](/icons/cvsweb/back.gif){kind=icon}
Return to face_enumerators_utilities.adb CVS log ![[TXT]](/icons/cvsweb/text.gif){kind=icon}
![[DIR]](/icons/cvsweb/dir.gif){kind=icon}
Up to [local] / OpenXM_contrib / PHC / Ada / Math_Lib / Supports