?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
fixed a small bug in eval ehrhart
[polylib.git] / doc / codeDoc / html / alpha_8c-source.html
blob8f23c600c15f4915054ec83a421fcbf3bf1cc0e4
1 <!DOCTYPE HTML PUBLIC "-//W3C//DTD HTML 4.0 Transitional//EN">
2 <html><head><meta http-equiv="Content-Type" content="text/html;charset=iso-8859-1">
3 <title>alpha.c Source File</title>
4 <link href="doxygen.css" rel="stylesheet" type="text/css">
5 </head><body>
6 <!-- Generated by Doxygen 1.2.15 -->
7 <center>
8 <a class="qindex" href="main.html">Main Page</a> &nbsp; <a class="qindex" href="annotated.html">Compound List</a> &nbsp; <a class="qindex" href="files.html">File List</a> &nbsp; <a class="qindex" href="functions.html">Compound Members</a> &nbsp; <a class="qindex" href="globals.html">File Members</a> &nbsp; </center>
9 <hr><h1>alpha.c</h1><a href="alpha_8c.html">Go to the documentation of this file.</a><div class="fragment"><pre>00001 <font class="comment">/* polytest.c */</font>
10 00002 <font class="preprocessor">#include &lt;stdio.h&gt;</font>
11 00003 <font class="preprocessor">#include &lt;polylib/polylib.h&gt;</font>
12 00004
13 00005 <font class="comment">/* alpha.c</font>
14 00006 <font class="comment"> COPYRIGHT</font>
15 00007 <font class="comment"> Both this software and its documentation are</font>
16 00008 <font class="comment"></font>
17 00009 <font class="comment"> Copyright 1993 by IRISA /Universite de Rennes I - France,</font>
18 00010 <font class="comment"> Copyright 1995,1996 by BYU, Provo, Utah</font>
19 00011 <font class="comment"> all rights reserved.</font>
20 00012 <font class="comment"></font>
21 00013 <font class="comment"> Permission is granted to copy, use, and distribute</font>
22 00014 <font class="comment"> for any commercial or noncommercial purpose under the terms</font>
23 00015 <font class="comment"> of the GNU General Public license, version 2, June 1991</font>
24 00016 <font class="comment"> (see file : LICENSING).</font>
25 00017 <font class="comment">*/</font>
26 00018
27 00019 <font class="preprocessor">#include &lt;stdio.h&gt;</font>
28 00020 <font class="preprocessor">#include &lt;stdlib.h&gt;</font>
29 00021 <font class="preprocessor">#include &lt;string.h&gt;</font>
30 00022
31 00023 <font class="comment">/*---------------------------------------------------------------------*/</font>
32 00024 <font class="comment">/* int exist_points(pos,P,context) */</font>
33 00025 <font class="comment">/* pos : index position of current loop index (0..hdim-1) */</font>
34 00026 <font class="comment">/* P: loop domain */</font>
35 00027 <font class="comment">/* context : context values for fixed indices */</font>
36 00028 <font class="comment">/* recursive procedure, recurs for each imbriquation */</font>
37 00029 <font class="comment">/* returns 1 if there exists any integer points in this polyhedron */</font>
38 00030 <font class="comment">/* returns 0 if there are none */</font>
39 00031 <font class="comment">/*---------------------------------------------------------------------*/</font>
40 <a name="l00032"></a><a class="code" href="alpha_8c.html#a0">00032</a> <font class="keyword">static</font> <font class="keywordtype">int</font> <a class="code" href="alpha_8c.html#a0">exist_points</a>(<font class="keywordtype">int</font> pos,Polyhedron *Pol,<a class="code" href="arithmetique_8h.html#a93">Value</a> *context) {
41 00033
42 00034 <a class="code" href="arithmetique_8h.html#a93">Value</a> LB, UB, k,tmp;
43 00035
44 00036 <a class="code" href="arithmetique_8h.html#a10">value_init</a>(LB); <a class="code" href="arithmetique_8h.html#a10">value_init</a>(UB);
45 00037 <a class="code" href="arithmetique_8h.html#a10">value_init</a>(k); <a class="code" href="arithmetique_8h.html#a10">value_init</a>(tmp);
46 00038 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(LB,0);
47 00039 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(UB,0);
48 00040
49 00041 <font class="comment">/* Problem if UB or LB is INFINITY */</font>
50 00042 <font class="keywordflow">if</font> (<a class="code" href="polyhedron_8c.html#a49">lower_upper_bounds</a>(pos,Pol,context,&amp;LB,&amp;UB) !=0) {
51 00043 <a class="code" href="errormsg_8c.html#a1">errormsg1</a>(<font class="stringliteral">"exist_points"</font>, <font class="stringliteral">"infdom"</font>, <font class="stringliteral">"infinite domain"</font>);
52 00044 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(LB);
53 00045 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(UB);
54 00046 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(k);
55 00047 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(tmp);
56 00048 <font class="keywordflow">return</font> -1;
57 00049 }
58 00050 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(context[pos],0);
59 00051 <font class="keywordflow">if</font>(<a class="code" href="arithmetique_8h.html#a26">value_lt</a>(UB,LB)) {
60 00052 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(LB);
61 00053 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(UB);
62 00054 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(k);
63 00055 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(tmp);
64 00056 <font class="keywordflow">return</font> 0;
65 00057 }
66 00058 <font class="keywordflow">if</font> (!Pol-&gt;next) {
67 00059 <a class="code" href="arithmetique_8h.html#a48">value_substract</a>(tmp,UB,LB);
68 00060 <a class="code" href="arithmetique_8h.html#a45">value_increment</a>(tmp,tmp);
69 00061 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(UB);
70 00062 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(LB);
71 00063 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(k);
72 00064 <font class="keywordflow">return</font> (<a class="code" href="arithmetique_8h.html#a63">value_pos_p</a>(tmp));
73 00065 }
74 00066
75 00067 <font class="keywordflow">for</font> (<a class="code" href="arithmetique_8h.html#a11">value_assign</a>(k,LB);<a class="code" href="arithmetique_8h.html#a27">value_le</a>(k,UB);<a class="code" href="arithmetique_8h.html#a45">value_increment</a>(k,k)) {
76 00068
77 00069 <font class="comment">/* insert k in context */</font>
78 00070 <a class="code" href="arithmetique_8h.html#a11">value_assign</a>(context[pos],k);
79 00071 <font class="keywordflow">if</font> (<a class="code" href="alpha_8c.html#a0">exist_points</a>(pos+1,Pol-&gt;next,context) &gt; 0 ) {
80 00072 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(LB); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(UB);
81 00073 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(k); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(tmp);
82 00074 <font class="keywordflow">return</font> 1;
83 00075 }
84 00076 }
85 00077 <font class="comment">/* Reset context */</font>
86 00078 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(context[pos],0);
87 00079 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(UB); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(LB);
88 00080 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(k); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(tmp);
89 00081 <font class="keywordflow">return</font> 0;
90 00082 }
91 00083
92 00084 <font class="comment">/*--------------------------------------------------------------*/</font>
93 00085 <font class="comment">/* Test to see if there are any integral points in a polyhedron */</font>
94 00086 <font class="comment">/*--------------------------------------------------------------*/</font>
95 <a name="l00087"></a><a class="code" href="alpha_8c.html#a1">00087</a> <font class="keywordtype">int</font> <a class="code" href="alpha_8c.html#a1">Polyhedron_Not_Empty</a>(Polyhedron *P,Polyhedron *C,<font class="keywordtype">int</font> <a class="code" href="verif__ehrhart_8c.html#a0">MAXRAYS</a>) {
96 00088
97 00089 <font class="keywordtype">int</font> res,i;
98 00090 <a class="code" href="arithmetique_8h.html#a93">Value</a> *context;
99 00091 Polyhedron *L;
100 00092
101 00093 <font class="comment">/* Create a context vector size dim+2 and set it to all zeros */</font>
102 00094 context = (<a class="code" href="arithmetique_8h.html#a93">Value</a> *) malloc((P-&gt;Dimension+2)*<font class="keyword">sizeof</font>(<a class="code" href="arithmetique_8h.html#a93">Value</a>));
103 00095
104 00096 <font class="comment">/* Initialize array 'context' */</font>
105 00097 <font class="keywordflow">for</font> (i=0;i&lt;(P-&gt;Dimension+2);i++)
106 00098 <a class="code" href="arithmetique_8h.html#a10">value_init</a>(context[i]);
107 00099
108 00100 <a class="code" href="vector_8c.html#a9">Vector_Set</a>(context,0,(P-&gt;Dimension+2));
109 00101
110 00102 <font class="comment">/* Set context[P-&gt;Dimension+1] = 1 (the constant) */</font>
111 00103 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(context[P-&gt;Dimension+1],1);
112 00104
113 00105 L = <a class="code" href="polyhedron_8c.html#a48">Polyhedron_Scan</a>(P,C,<a class="code" href="verif__ehrhart_8c.html#a0">MAXRAYS</a>);
114 00106 res = <a class="code" href="alpha_8c.html#a0">exist_points</a>(1,L,context);
115 00107 <a class="code" href="polyhedron_8c.html#a19">Domain_Free</a>(L);
116 00108
117 00109 <font class="comment">/* Clear array 'context' */</font>
118 00110 <font class="keywordflow">for</font> (i=0;i&lt;(P-&gt;Dimension+2);i++)
119 00111 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(context[i]);
120 00112 free(context);
121 00113 <font class="keywordflow">return</font> res;
122 00114 }
123 00115
124 00116 <font class="comment">/* PolyhedronLTQ ( P1, P2 ) */</font>
125 00117 <font class="comment">/* P1 and P2 must be simple polyhedra */</font>
126 00118 <font class="comment">/* result = 0 --&gt; not comparable */</font>
127 00119 <font class="comment">/* result = -1 --&gt; P1 &lt; P2 */</font>
128 00120 <font class="comment">/* result = 1 --&gt; P1 &gt; P2 */</font>
129 00121 <font class="comment">/* INDEX = 1 .... Dimension */</font>
130 <a name="l00122"></a><a class="code" href="alpha_8c.html#a2">00122</a> <font class="keywordtype">int</font> <a class="code" href="alpha_8c.html#a2">PolyhedronLTQ</a> (Polyhedron *Pol1,Polyhedron *Pol2,<font class="keywordtype">int</font> INDEX, <font class="keywordtype">int</font> PDIM, <font class="keywordtype">int</font> NbMaxConstrs) {
131 00123
132 00124 <font class="keywordtype">int</font> res, dim, i, j, k;
133 00125 Polyhedron *Q1, *Q2, *Q3, *Q4, *Q;
134 00126 Matrix *Mat;
135 00127
136 00128 <font class="keywordflow">if</font> (Pol1-&gt;next || Pol2-&gt;next) {
137 00129 <a class="code" href="errormsg_8c.html#a1">errormsg1</a>(<font class="stringliteral">"PolyhedronLTQ"</font>, <font class="stringliteral">"compoly"</font>, <font class="stringliteral">"Can only compare polyhedra"</font>);
138 00130 <font class="keywordflow">return</font> 0;
139 00131 }
140 00132 <font class="keywordflow">if</font> (Pol1-&gt;Dimension != Pol2-&gt;Dimension) {
141 00133 <a class="code" href="errormsg_8c.html#a1">errormsg1</a>(<font class="stringliteral">"PolyhedronLTQ"</font>,<font class="stringliteral">"diffdim"</font>,<font class="stringliteral">"Polyhedra are not same dimension"</font>);
142 00134 <font class="keywordflow">return</font> 0;
143 00135 }
144 00136 dim = Pol1-&gt;Dimension+2;
145 00137
146 00138 <font class="preprocessor">#ifdef DEBUG</font>
147 00139 <font class="preprocessor"></font> fprintf(stdout, <font class="stringliteral">"P1\n"</font>);
148 00140 <a class="code" href="polyhedron_8c.html#a20">Polyhedron_Print</a>(stdout,P_VALUE_FMT,Pol1);
149 00141 fprintf(stdout, <font class="stringliteral">"P2\n"</font>);
150 00142 <a class="code" href="polyhedron_8c.html#a20">Polyhedron_Print</a>(stdout,P_VALUE_FMT,Pol2);
151 00143 <font class="preprocessor">#endif</font>
152 00144 <font class="preprocessor"></font>
153 00145 <font class="comment">/* Create the Line to add */</font>
154 00146 k = Pol1-&gt;Dimension-INDEX+1-PDIM;
155 00147 Mat = <a class="code" href="matrix_8c.html#a0">Matrix_Alloc</a>(k,dim);
156 00148 <a class="code" href="vector_8c.html#a9">Vector_Set</a>(Mat-&gt;p_Init,0,dim*k);
157 00149 <font class="keywordflow">for</font>(j=0,i=INDEX;j&lt;k;i++,j++)
158 00150 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(Mat-&gt;p[j][i],1);
159 00151
160 00152 Q1 = <a class="code" href="polyhedron_8c.html#a34">AddRays</a>(Mat-&gt;p[0],k,Pol1,NbMaxConstrs);
161 00153 Q2 = <a class="code" href="polyhedron_8c.html#a34">AddRays</a>(Mat-&gt;p[0],k,Pol2,NbMaxConstrs);
162 00154
163 00155 <font class="preprocessor">#ifdef DEBUG</font>
164 00156 <font class="preprocessor"></font> fprintf(stdout, <font class="stringliteral">"Q1\n"</font>);
165 00157 <a class="code" href="polyhedron_8c.html#a20">Polyhedron_Print</a>(stdout,P_VALUE_FMT,Q1);
166 00158 fprintf(stdout, <font class="stringliteral">"Q2\n"</font>);
167 00159 <a class="code" href="polyhedron_8c.html#a20">Polyhedron_Print</a>(stdout,P_VALUE_FMT,Q2);
168 00160 <font class="preprocessor">#endif</font>
169 00161 <font class="preprocessor"></font>
170 00162 <a class="code" href="matrix_8c.html#a1">Matrix_Free</a>(Mat);
171 00163 Q = <a class="code" href="polyhedron_8c.html#a32">DomainIntersection</a>(Q1,Q2,NbMaxConstrs);
172 00164
173 00165 <font class="preprocessor">#ifdef DEBUG</font>
174 00166 <font class="preprocessor"></font> fprintf(stdout, <font class="stringliteral">"Q\n"</font>);
175 00167 <a class="code" href="polyhedron_8c.html#a20">Polyhedron_Print</a>(stdout,P_VALUE_FMT,Q);
176 00168 <font class="preprocessor">#endif</font>
177 00169 <font class="preprocessor"></font>
178 00170 <a class="code" href="polyhedron_8c.html#a19">Domain_Free</a>(Q1);
179 00171 <a class="code" href="polyhedron_8c.html#a19">Domain_Free</a>(Q2);
180 00172
181 00173 <font class="keywordflow">if</font> (emptyQ(Q)) res = 0; <font class="comment">/* not comparable */</font>
182 00174 <font class="keywordflow">else</font> {
183 00175 Q1 = <a class="code" href="polyhedron_8c.html#a32">DomainIntersection</a>(Pol1,Q,NbMaxConstrs);
184 00176 Q2 = <a class="code" href="polyhedron_8c.html#a32">DomainIntersection</a>(Pol2,Q,NbMaxConstrs);
185 00177
186 00178 <font class="preprocessor">#ifdef DEBUG</font>
187 00179 <font class="preprocessor"></font> fprintf(stdout, <font class="stringliteral">"Q1\n"</font>);
188 00180 <a class="code" href="polyhedron_8c.html#a20">Polyhedron_Print</a>(stdout,P_VALUE_FMT,Q1);
189 00181 fprintf(stdout, <font class="stringliteral">"Q2\n"</font>);
190 00182 <a class="code" href="polyhedron_8c.html#a20">Polyhedron_Print</a>(stdout,P_VALUE_FMT,Q2);
191 00183 <font class="preprocessor">#endif</font>
192 00184 <font class="preprocessor"></font>
193 00185 k = Q1-&gt;NbConstraints + Q2-&gt;NbConstraints;
194 00186 Mat = <a class="code" href="matrix_8c.html#a0">Matrix_Alloc</a>(k, dim);
195 00187 <a class="code" href="vector_8c.html#a9">Vector_Set</a>(Mat-&gt;p_Init,0,k*dim);
196 00188
197 00189 <font class="comment">/* First compute surrounding polyhedron */</font>
198 00190 j=0;
199 00191 <font class="keywordflow">for</font> (i=0; i&lt;Q1-&gt;NbConstraints; i++) {
200 00192 <font class="keywordflow">if</font> ((<a class="code" href="arithmetique_8h.html#a69">value_one_p</a>(Q1-&gt;Constraint[i][0])) &amp;&amp; (<a class="code" href="arithmetique_8h.html#a63">value_pos_p</a>(Q1-&gt;Constraint[i][INDEX]))) {
201 00193
202 00194 <font class="comment">/* keep Q1's lower bounds */</font>
203 00195 <font class="keywordflow">for</font> (k=0; k&lt;dim; k++)
204 00196 <a class="code" href="arithmetique_8h.html#a11">value_assign</a>(Mat-&gt;p[j][k],Q1-&gt;Constraint[i][k]);
205 00197 j++;
206 00198 }
207 00199 }
208 00200 <font class="keywordflow">for</font> (i=0; i&lt;Q2-&gt;NbConstraints; i++) {
209 00201 <font class="keywordflow">if</font> ((<a class="code" href="arithmetique_8h.html#a69">value_one_p</a>(Q2-&gt;Constraint[i][0])) &amp;&amp; (<a class="code" href="arithmetique_8h.html#a64">value_neg_p</a>(Q2-&gt;Constraint[i][INDEX]))) {
210 00202
211 00203 <font class="comment">/* and keep Q2's upper bounds */</font>
212 00204 <font class="keywordflow">for</font> (k=0; k&lt;dim; k++)
213 00205 <a class="code" href="arithmetique_8h.html#a11">value_assign</a>(Mat-&gt;p[j][k],Q2-&gt;Constraint[i][k]);
214 00206 j++;
215 00207 }
216 00208 }
217 00209 Q4 = <a class="code" href="polyhedron_8c.html#a28">AddConstraints</a>(Mat-&gt;p[0], j, Q, NbMaxConstrs);
218 00210 <a class="code" href="matrix_8c.html#a1">Matrix_Free</a>(Mat);
219 00211
220 00212 <font class="preprocessor">#ifdef debug</font>
221 00213 <font class="preprocessor"></font> fprintf(stderr, <font class="stringliteral">"Q4 surrounding polyhedron\n"</font>);
222 00214 Polyhderon_Print(stderr,P_VALUE_FMT, Q4);
223 00215 <font class="preprocessor">#endif</font>
224 00216 <font class="preprocessor"></font>
225 00217 <font class="comment">/* if surrounding polyhedron is empty, D1&gt;D2 */</font>
226 00218 <font class="keywordflow">if</font> (emptyQ(Q4)) {
227 00219 res = 1;
228 00220
229 00221 <font class="preprocessor">#ifdef debug</font>
230 00222 <font class="preprocessor"></font> fprintf(stderr, <font class="stringliteral">"Surrounding polyhedron is empty\n"</font>);
231 00223 <font class="preprocessor">#endif</font>
232 00224 <font class="preprocessor"></font> <font class="keywordflow">goto</font> LTQdone2;
233 00225 }
234 00226
235 00227 <font class="comment">/* Test if Q1 &lt; Q2 */</font>
236 00228 <font class="comment">/* Build a constraint array for &gt;= Q1 and &lt;= Q2 */</font>
237 00229 Mat = <a class="code" href="matrix_8c.html#a0">Matrix_Alloc</a>(2,dim);
238 00230 <a class="code" href="vector_8c.html#a9">Vector_Set</a>(Mat-&gt;p_Init,0,2*dim);
239 00231
240 00232 <font class="comment">/* Choose a contraint from Q1 */</font>
241 00233 <font class="keywordflow">for</font> (i=0; i&lt;Q1-&gt;NbConstraints; i++) {
242 00234 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a67">value_zero_p</a>(Q1-&gt;Constraint[i][0])) {
243 00235
244 00236 <font class="comment">/* Equality */</font>
245 00237 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a67">value_zero_p</a>(Q1-&gt;Constraint[i][INDEX])) {
246 00238
247 00239 <font class="comment">/* Ignore side constraint (they are in Q) */</font>
248 00240 <font class="keywordflow">continue</font>;
249 00241 }
250 00242 <font class="keywordflow">else</font> <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a64">value_neg_p</a>(Q1-&gt;Constraint[i][INDEX])) {
251 00243
252 00244 <font class="comment">/* copy -constraint to Mat */</font>
253 00245 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(Mat-&gt;p[0][0],1);
254 00246 <font class="keywordflow">for</font> (k=1; k&lt;dim; k++)
255 00247 <a class="code" href="arithmetique_8h.html#a54">value_oppose</a>(Mat-&gt;p[0][k],Q1-&gt;Constraint[i][k]);
256 00248 }
257 00249 <font class="keywordflow">else</font> {
258 00250
259 00251 <font class="comment">/* Copy constraint to Mat */</font>
260 00252
261 00253 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(Mat-&gt;p[0][0],1);
262 00254 <font class="keywordflow">for</font> (k=1; k&lt;dim; k++)
263 00255 <a class="code" href="arithmetique_8h.html#a11">value_assign</a>(Mat-&gt;p[0][k],Q1-&gt;Constraint[i][k]);
264 00256 }
265 00257 }
266 00258 <font class="keywordflow">else</font> <font class="keywordflow">if</font>(<a class="code" href="arithmetique_8h.html#a64">value_neg_p</a>(Q1-&gt;Constraint[i][INDEX])) {
267 00259
268 00260 <font class="comment">/* Upper bound -- make a lower bound from it */</font>
269 00261 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(Mat-&gt;p[0][0],1);
270 00262 <font class="keywordflow">for</font> (k=1; k&lt;dim; k++)
271 00263 <a class="code" href="arithmetique_8h.html#a54">value_oppose</a>(Mat-&gt;p[0][k],Q1-&gt;Constraint[i][k]);
272 00264 }
273 00265 <font class="keywordflow">else</font> {
274 00266
275 00267 <font class="comment">/* Lower or side bound -- ignore it */</font>
276 00268 <font class="keywordflow">continue</font>;
277 00269 }
278 00270
279 00271 <font class="comment">/* Choose a constraint from Q2 */</font>
280 00272 <font class="keywordflow">for</font> (j=0; j&lt;Q2-&gt;NbConstraints; j++) {
281 00273 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a67">value_zero_p</a>(Q2-&gt;Constraint[j][0])) { <font class="comment">/* equality */</font>
282 00274 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a67">value_zero_p</a>(Q2-&gt;Constraint[j][INDEX])) {
283 00275
284 00276 <font class="comment">/* Ignore side constraint (they are in Q) */</font>
285 00277 <font class="keywordflow">continue</font>;
286 00278 }
287 00279 <font class="keywordflow">else</font> <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a63">value_pos_p</a>(Q2-&gt;Constraint[j][INDEX])) {
288 00280
289 00281 <font class="comment">/* Copy -constraint to Mat */</font>
290 00282 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(Mat-&gt;p[1][0],1);
291 00283 <font class="keywordflow">for</font> (k=1; k&lt;dim; k++)
292 00284 <a class="code" href="arithmetique_8h.html#a54">value_oppose</a>(Mat-&gt;p[1][k],Q2-&gt;Constraint[j][k]);
293 00285 }
294 00286 <font class="keywordflow">else</font> {
295 00287
296 00288 <font class="comment">/* Copy constraint to Mat */</font>
297 00289 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(Mat-&gt;p[1][0],1);
298 00290 <font class="keywordflow">for</font> (k=1; k&lt;dim; k++)
299 00291 <a class="code" href="arithmetique_8h.html#a11">value_assign</a>(Mat-&gt;p[1][k],Q2-&gt;Constraint[j][k]);
300 00292 };
301 00293 }
302 00294 <font class="keywordflow">else</font> <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a63">value_pos_p</a>(Q2-&gt;Constraint[j][INDEX])) {
303 00295
304 00296 <font class="comment">/* Lower bound -- make an upper bound from it */</font>
305 00297 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(Mat-&gt;p[1][0],1);
306 00298 <font class="keywordflow">for</font>(k=1;k&lt;dim;k++)
307 00299 <a class="code" href="arithmetique_8h.html#a54">value_oppose</a>(Mat-&gt;p[1][k],Q2-&gt;Constraint[j][k]);
308 00300 }
309 00301 <font class="keywordflow">else</font> {
310 00302
311 00303 <font class="comment">/* Upper or side bound -- ignore it */</font>
312 00304 <font class="keywordflow">continue</font>;
313 00305 };
314 00306
315 00307 <font class="preprocessor">#ifdef DEBUG</font>
316 00308 <font class="preprocessor"></font> fprintf(stdout, <font class="stringliteral">"i=%d j=%d M=\n"</font>, i+1, j+1);
317 00309 <a class="code" href="matrix_8c.html#a2">Matrix_Print</a>(stdout,P_VALUE_FMT,Mat);
318 00310 <font class="preprocessor">#endif</font>
319 00311 <font class="preprocessor"></font>
320 00312 <font class="comment">/* Add Mat to Q and see if anything is made */</font>
321 00313 Q3 = <a class="code" href="polyhedron_8c.html#a28">AddConstraints</a>(Mat-&gt;p[0],2,Q,NbMaxConstrs);
322 00314
323 00315 <font class="preprocessor">#ifdef DEBUG</font>
324 00316 <font class="preprocessor"></font> fprintf(stdout, <font class="stringliteral">"Q3\n"</font>);
325 00317 <a class="code" href="polyhedron_8c.html#a20">Polyhedron_Print</a>(stdout,P_VALUE_FMT,Q3);
326 00318 <font class="preprocessor">#endif</font>
327 00319 <font class="preprocessor"></font>
328 00320 <font class="keywordflow">if</font> (!emptyQ(Q3)) {
329 00321 <a class="code" href="polyhedron_8c.html#a19">Domain_Free</a>(Q3);
330 00322
331 00323 <font class="preprocessor">#ifdef DEBUG</font>
332 00324 <font class="preprocessor"></font> fprintf(stdout, <font class="stringliteral">"not empty\n"</font>);
333 00325 <font class="preprocessor">#endif</font>
334 00326 <font class="preprocessor"></font> res = -1;
335 00327 <font class="keywordflow">goto</font> LTQdone;
336 00328 }
337 00329 <font class="preprocessor">#ifdef DEBUG</font>
338 00330 <font class="preprocessor"></font> fprintf(stdout,<font class="stringliteral">"empty\n"</font>);
339 00331 <font class="preprocessor">#endif</font>
340 00332 <font class="preprocessor"></font> <a class="code" href="polyhedron_8c.html#a19">Domain_Free</a>(Q3);
341 00333 } <font class="comment">/* end for j */</font>
342 00334 } <font class="comment">/* end for i */</font>
343 00335 res = 1;
344 00336 LTQdone:
345 00337 <a class="code" href="matrix_8c.html#a1">Matrix_Free</a>(Mat);
346 00338 LTQdone2:
347 00339 <a class="code" href="polyhedron_8c.html#a19">Domain_Free</a>(Q4);
348 00340 <a class="code" href="polyhedron_8c.html#a19">Domain_Free</a>(Q1);
349 00341 <a class="code" href="polyhedron_8c.html#a19">Domain_Free</a>(Q2);
350 00342 }
351 00343 <a class="code" href="polyhedron_8c.html#a19">Domain_Free</a>(Q);
352 00344
353 00345 <font class="preprocessor">#ifdef DEBUG</font>
354 00346 <font class="preprocessor"></font> fprintf(stdout, <font class="stringliteral">"res = %d\n"</font>, res);
355 00347 <font class="preprocessor">#endif</font>
356 00348 <font class="preprocessor"></font>
357 00349 <font class="keywordflow">return</font> res;
358 00350 } <font class="comment">/* PolyhedronLTQ */</font>
359 00351
360 00352 <font class="comment">/* GaussSimplify --</font>
361 00353 <font class="comment"> Given Mat1, a matrix of equalities, performs Gaussian elimination.</font>
362 00354 <font class="comment"> Find a minimum basis, Returns the rank.</font>
363 00355 <font class="comment"> Mat1 is context, Mat2 is reduced in context of Mat1</font>
364 00356 <font class="comment">*/</font>
365 <a name="l00357"></a><a class="code" href="alpha_8c.html#a3">00357</a> <font class="keywordtype">int</font> <a class="code" href="alpha_8c.html#a3">GaussSimplify</a>(Matrix *Mat1,Matrix *Mat2) {
366 00358
367 00359 <font class="keywordtype">int</font> NbRows = Mat1-&gt;NbRows;
368 00360 <font class="keywordtype">int</font> NbCols = Mat1-&gt;NbColumns;
369 00361 <font class="keywordtype">int</font> *column_index;
370 00362 <font class="keywordtype">int</font> i, j, k, n, t, pivot, Rank;
371 00363 <a class="code" href="arithmetique_8h.html#a93">Value</a> gcd, tmp, *cp;
372 00364
373 00365 column_index=(<font class="keywordtype">int</font> *)malloc(NbCols * <font class="keyword">sizeof</font>(<font class="keywordtype">int</font>));
374 00366 <font class="keywordflow">if</font> (!column_index) {
375 00367 <a class="code" href="errormsg_8c.html#a1">errormsg1</a>(<font class="stringliteral">"GaussSimplify"</font>, <font class="stringliteral">"outofmem"</font>, <font class="stringliteral">"out of memory space\n"</font>);
376 00368 Pol_status = 1;
377 00369 <font class="keywordflow">return</font> 0;
378 00370 }
379 00371
380 00372 <font class="comment">/* Initialize all the 'Value' variables */</font>
381 00373 <a class="code" href="arithmetique_8h.html#a10">value_init</a>(gcd); <a class="code" href="arithmetique_8h.html#a10">value_init</a>(tmp);
382 00374
383 00375 Rank=0;
384 00376 <font class="keywordflow">for</font> (j=0; j&lt;NbCols; j++) { <font class="comment">/* for each column starting at */</font>
385 00377 <font class="keywordflow">for</font> (i=Rank; i&lt;NbRows; i++) <font class="comment">/* diagonal, look down to find */</font>
386 00378 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a68">value_notzero_p</a>(Mat1-&gt;p[i][j])) <font class="comment">/* the first non-zero entry */</font>
387 00379 <font class="keywordflow">break</font>;
388 00380 <font class="keywordflow">if</font> (i!=NbRows) { <font class="comment">/* was one found ? */</font>
389 00381 <font class="keywordflow">if</font> (i!=Rank) <font class="comment">/* was it found below the diagonal?*/</font>
390 00382 <a class="code" href="vector_8c.html#a10">Vector_Exchange</a>(Mat1-&gt;p[Rank],Mat1-&gt;p[i],NbCols);
391 00383
392 00384 <font class="comment">/* Normalize the pivot row */</font>
393 00385 <a class="code" href="vector_8c.html#a23">Vector_Gcd</a>(Mat1-&gt;p[Rank],NbCols,&amp;gcd);
394 00386
395 00387 <font class="comment">/* If (gcd &gt;= 2) */</font>
396 00388 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(tmp,2);
397 00389 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a25">value_ge</a>(gcd,tmp)) {
398 00390 cp = Mat1-&gt;p[Rank];
399 00391 <font class="keywordflow">for</font> (k=0; k&lt;NbCols; k++,cp++)
400 00392 <a class="code" href="arithmetique_8h.html#a51">value_division</a>(*cp,*cp,gcd);
401 00393 }
402 00394 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a64">value_neg_p</a>(Mat1-&gt;p[Rank][j])) {
403 00395 cp = Mat1-&gt;p[Rank];
404 00396 <font class="keywordflow">for</font> (k=0; k&lt;NbCols; k++,cp++)
405 00397 <a class="code" href="arithmetique_8h.html#a54">value_oppose</a>(*cp,*cp);
406 00398 }
407 00399 <font class="comment">/* End of normalize */</font>
408 00400 pivot=i;
409 00401 <font class="keywordflow">for</font> (i=0;i&lt;NbRows;i++) <font class="comment">/* Zero out the rest of the column */</font>
410 00402 <font class="keywordflow">if</font> (i!=Rank) {
411 00403 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a68">value_notzero_p</a>(Mat1-&gt;p[i][j])) {
412 00404 <a class="code" href="arithmetique_8h.html#a93">Value</a> a, a1, a2, a1abs, a2abs;
413 00405 <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a); <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a1); <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a2);
414 00406 <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a1abs); <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a2abs);
415 00407 <a class="code" href="arithmetique_8h.html#a11">value_assign</a>(a1,Mat1-&gt;p[i][j]);
416 00408 <a class="code" href="arithmetique_8h.html#a55">value_absolute</a>(a1abs,a1);
417 00409 <a class="code" href="arithmetique_8h.html#a11">value_assign</a>(a2,Mat1-&gt;p[Rank][j]);
418 00410 <a class="code" href="arithmetique_8h.html#a55">value_absolute</a>(a2abs,a2);
419 00411 <a class="code" href="vector_8c.html#a3">Gcd</a>(a1abs,a2abs,&amp;a);
420 00412 <a class="code" href="arithmetique_8h.html#a51">value_division</a>(a1,a1,a);
421 00413 <a class="code" href="arithmetique_8h.html#a51">value_division</a>(a2,a2,a);
422 00414 <a class="code" href="arithmetique_8h.html#a54">value_oppose</a>(a1,a1);
423 00415 <a class="code" href="vector_8c.html#a20">Vector_Combine</a>(Mat1-&gt;p[i],Mat1-&gt;p[Rank],Mat1-&gt;p[i],a2,
424 00416 a1,NbCols);
425 00417 <a class="code" href="vector_8c.html#a25">Vector_Normalize</a>(Mat1-&gt;p[i],NbCols);
426 00418 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a1); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a2);
427 00419 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a1abs); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a2abs);
428 00420 }
429 00421 }
430 00422 column_index[Rank]=j;
431 00423 Rank++;
432 00424 }
433 00425 } <font class="comment">/* end of Gauss elimination */</font>
434 00426
435 00427
436 00428 <font class="keywordflow">if</font> (Mat2) { <font class="comment">/* Mat2 is a transformation matrix (i,j-&gt;f(i,j))....</font>
437 00429 <font class="comment"> can't scale it because can't scale both sides of -&gt; */</font>
438 00430 <font class="comment">/* normalizes an affine transformation */</font>
439 00431 <font class="comment">/* priority of forms */</font>
440 00432 <font class="comment">/* 1. i' -&gt; i (identity) */</font>
441 00433 <font class="comment">/* 2. i' -&gt; i + constant (uniform) */</font>
442 00434 <font class="comment">/* 3. i' -&gt; constant (broadcast) */</font>
443 00435 <font class="comment">/* 4. i' -&gt; j (permutation) */</font>
444 00436 <font class="comment">/* 5. i' -&gt; j + constant ( ) */</font>
445 00437 <font class="comment">/* 6. i' -&gt; i + j + constant (non-uniform) */</font>
446 00438 <font class="keywordflow">for</font> (k=0; k&lt;Rank; k++) {
447 00439 j = column_index[k];
448 00440 <font class="keywordflow">for</font> (i=0; i&lt;(Mat2-&gt;NbRows-1);i++) { <font class="comment">/* all but the last row 0...0 1 */</font>
449 00441 <font class="keywordflow">if</font> ((i!=j) &amp;&amp; <a class="code" href="arithmetique_8h.html#a68">value_notzero_p</a>(Mat2-&gt;p[i][j])) {
450 00442
451 00443 <font class="comment">/* Remove dependency of i' on j */</font>
452 00444 <a class="code" href="arithmetique_8h.html#a93">Value</a> a, a1, a1abs, a2, a2abs;
453 00445 <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a); <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a1); <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a2);
454 00446 <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a1abs); <a class="code" href="arithmetique_8h.html#a10">value_init</a>(a2abs);
455 00447 <a class="code" href="arithmetique_8h.html#a11">value_assign</a>(a1,Mat2-&gt;p[i][j]);
456 00448 <a class="code" href="arithmetique_8h.html#a55">value_absolute</a>(a1abs,a1);
457 00449 <a class="code" href="arithmetique_8h.html#a11">value_assign</a>(a2,Mat1-&gt;p[k][j]);
458 00450 <a class="code" href="arithmetique_8h.html#a55">value_absolute</a>(a2abs,a2);
459 00451 <a class="code" href="vector_8c.html#a3">Gcd</a>(a1abs,a2abs,&amp;a);
460 00452 <a class="code" href="arithmetique_8h.html#a51">value_division</a>(a1,a1,a);
461 00453 <a class="code" href="arithmetique_8h.html#a51">value_division</a>(a2,a2,a);
462 00454 <a class="code" href="arithmetique_8h.html#a54">value_oppose</a>(a1,a1);
463 00455 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a69">value_one_p</a>(a2)) {
464 00456 <a class="code" href="vector_8c.html#a20">Vector_Combine</a>(Mat2-&gt;p[i],Mat1-&gt;p[k],Mat2-&gt;p[i],a2,
465 00457 a1,NbCols);
466 00458
467 00459 <font class="comment">/* Vector_Normalize(Mat2-&gt;p[i],NbCols); -- can't do T */</font>
468 00460 } <font class="comment">/* otherwise, can't do it without mult lhs prod (2i,3j-&gt;...) */</font>
469 00461 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a1); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a2);
470 00462 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a1abs); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(a2abs);
471 00463
472 00464 }
473 00465 <font class="keywordflow">else</font> <font class="keywordflow">if</font> ((i==j) &amp;&amp; <a class="code" href="arithmetique_8h.html#a67">value_zero_p</a>(Mat2-&gt;p[i][j])) {
474 00466
475 00467 <font class="comment">/* 'i' does not depend on j */</font>
476 00468 <font class="keywordflow">for</font> (n=j+1; n &lt; (NbCols-1); n++) {
477 00469 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a68">value_notzero_p</a>(Mat2-&gt;p[i][n])) { <font class="comment">/* i' depends on some n */</font>
478 00470 <a class="code" href="arithmetique_8h.html#a12">value_set_si</a>(tmp,1);
479 00471 <a class="code" href="vector_8c.html#a20">Vector_Combine</a>(Mat2-&gt;p[i],Mat1-&gt;p[k],Mat2-&gt;p[i],tmp,
480 00472 tmp,NbCols);
481 00473 <font class="keywordflow">break</font>;
482 00474 } <font class="comment">/* if 'i' depends on just a constant, then leave it alone.*/</font>
483 00475 }
484 00476 }
485 00477 }
486 00478 }
487 00479
488 00480 <font class="comment">/* Check last row of transformation Mat2 */</font>
489 00481 <font class="keywordflow">for</font> (j=0; j&lt;(NbCols-1); j++)
490 00482 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a68">value_notzero_p</a>(Mat2-&gt;p[Mat2-&gt;NbRows-1][j])) {
491 00483 <a class="code" href="errormsg_8c.html#a1">errormsg1</a>(<font class="stringliteral">"GaussSimplify"</font>, <font class="stringliteral">"corrtrans"</font>, <font class="stringliteral">"Corrupted transformation\n"</font>);
492 00484 <font class="keywordflow">break</font>;
493 00485 }
494 00486
495 00487 <font class="keywordflow">if</font> (<a class="code" href="arithmetique_8h.html#a70">value_notone_p</a>(Mat2-&gt;p[Mat2-&gt;NbRows-1][NbCols-1])) {
496 00488 <a class="code" href="errormsg_8c.html#a1">errormsg1</a>(<font class="stringliteral">"GaussSimplify"</font>, <font class="stringliteral">"corrtrans"</font>, <font class="stringliteral">"Corrupted transformation\n"</font>);
497 00489 }
498 00490 }
499 00491 <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(gcd); <a class="code" href="arithmetique_8h.html#a14">value_clear</a>(tmp);
500 00492 free(column_index);
501 00493 <font class="keywordflow">return</font> Rank;
502 00494 } <font class="comment">/* GaussSimplify */</font>
503 00495
504 00496 <font class="comment">/* </font>
505 00497 <font class="comment"> * Topologically sort 'n' polyhdera and return 0 on failure, otherwise return </font>
506 00498 <font class="comment"> * 1 on success. Here 'L' is a an array of pointers to polyhedra, 'n' is the </font>
507 00499 <font class="comment"> * number of polyhedra, 'index' is the level to consider for partial ordering</font>
508 00500 <font class="comment"> * 'pdim' is the parameter space dimension, 'time' is an array of 'n' integers</font>
509 00501 <font class="comment"> * to store logical time values, 'pvect', if not NULL, is an array of 'n' </font>
510 00502 <font class="comment"> * integers that contains a permutation specification after call and MAXRAYS is</font>
511 00503 <font class="comment"> * the workspace size for polyhedra operations. </font>
512 00504 <font class="comment"> */</font>
513 <a name="l00505"></a><a class="code" href="alpha_8c.html#a4">00505</a> <font class="keywordtype">int</font> <a class="code" href="alpha_8c.html#a4">PolyhedronTSort</a> (Polyhedron **L,<font class="keywordtype">unsigned</font> <font class="keywordtype">int</font> n,<font class="keywordtype">unsigned</font> <font class="keywordtype">int</font> index,<font class="keywordtype">unsigned</font> <font class="keywordtype">int</font> pdim,<font class="keywordtype">int</font> *time,<font class="keywordtype">int</font> *pvect,<font class="keywordtype">unsigned</font> <font class="keywordtype">int</font> <a class="code" href="verif__ehrhart_8c.html#a0">MAXRAYS</a>) {
514 00506
515 00507 <font class="keywordtype">unsigned</font> <font class="keywordtype">int</font> <font class="keyword">const</font> nbcells = ((n*(n-1))&gt;&gt;1); <font class="comment">/* Number of memory cells </font>
516 00508 <font class="comment"> to allocate, see below */</font>
517 00509 <font class="keywordtype">int</font> *dag; <font class="comment">/* The upper triangular matrix */</font>
518 00510 <font class="keywordtype">int</font> **p; <font class="comment">/* Array of matrix row addresses */</font>
519 00511 <font class="keywordtype">unsigned</font> <font class="keywordtype">int</font> i, j, k;
520 00512 <font class="keywordtype">unsigned</font> <font class="keywordtype">int</font> t, nb, isroot;
521 00513
522 00514 <font class="keywordflow">if</font> (n&lt;2) <font class="keywordflow">return</font> 0;
523 00515
524 00516 <font class="comment">/* we need an upper triangular matrix (example with n=4):</font>
525 00517 <font class="comment"></font>
526 00518 <font class="comment"> . o o o</font>
527 00519 <font class="comment"> . . o o . are unuseful cells, o are useful cells</font>
528 00520 <font class="comment"> . . . o</font>
529 00521 <font class="comment"> . . . .</font>
530 00522 <font class="comment"> </font>
531 00523 <font class="comment"> so we need to allocate (n)(n-1)/2 cells</font>
532 00524 <font class="comment"> - dag will point to this memory.</font>
533 00525 <font class="comment"> - p[i] will point to row i of the matrix</font>
534 00526 <font class="comment"> p[0] = dag - 1 (such that p[0][1] == dag[0])</font>
535 00527 <font class="comment"> p[1] = dag - 1 + (n-1)</font>
536 00528 <font class="comment"> p[2] = dag - 1 + (n-1) + (n-2)</font>
537 00529 <font class="comment"> ...</font>
538 00530 <font class="comment"> p[i] = p[i-1] + (n-1-i)</font>
539 00531 <font class="comment"> */</font>
540 00532
541 00533 <font class="comment">/* malloc n(n-1)/2 for dag and n-1 for p (node n does not have any row) */</font>
542 00534 dag = (<font class="keywordtype">int</font> *) malloc(nbcells*<font class="keyword">sizeof</font>(<font class="keywordtype">int</font>));
543 00535 <font class="keywordflow">if</font> (!dag) <font class="keywordflow">return</font> 0;
544 00536 p = (<font class="keywordtype">int</font> **) malloc ((n-1) * <font class="keyword">sizeof</font>(<font class="keywordtype">int</font> *));
545 00537 <font class="keywordflow">if</font> (!p) {
546 00538 free(dag); <font class="keywordflow">return</font> 0;
547 00539 }
548 00540
549 00541 <font class="comment">/* Initialize 'p' and 'dag' */</font>
550 00542 p[0] = dag-1;
551 00543 <font class="keywordflow">for</font> (i=1; i&lt;n-1; i++)
552 00544 p[i] = p[i-1] + (n-1)-i;
553 00545 <font class="keywordflow">for</font> (i=0; i&lt;nbcells; i++)
554 00546 dag[i] = -2; <font class="comment">/* -2 means 'not computed yet' */</font>
555 00547 <font class="keywordflow">for</font> (i=0; i&lt;n; i++) time[i] = -1;
556 00548
557 00549 <font class="comment">/* Compute the dag using transitivity to reduce the number of */</font>
558 00550 <font class="comment">/* PolyhedronLTQ calls. */</font>
559 00551 <font class="keywordflow">for</font> (i=0; i&lt;n-1; i++) {
560 00552 <font class="keywordflow">for</font> (j=i+1; j&lt;n; j++) {
561 00553 <font class="keywordflow">if</font> (p[i][j] == -2) <font class="comment">/* not computed yes */</font>
562 00554 p[i][j] = <a class="code" href="alpha_8c.html#a2">PolyhedronLTQ</a>(L[i], L[j], index, pdim, <a class="code" href="verif__ehrhart_8c.html#a0">MAXRAYS</a>);
563 00555 <font class="keywordflow">if</font> (p[i][j] != 0) {
564 00556
565 00557 <font class="comment">/* if p[i][j] is 1 or -1, look for -p[i][j] on the same row:</font>
566 00558 <font class="comment"> p[i][j] == -p[i][k] ==&gt; p[j][k] = p[i][k] (transitivity)</font>
567 00559 <font class="comment"> note: p[r][c] == -p[c][r], use this to avoid reading or writing</font>
568 00560 <font class="comment"> under the matrix diagonal</font>
569 00561 <font class="comment"> */</font>
570 00562
571 00563 <font class="comment">/* first, k&lt;i so look for p[i][j] == p[k][i] (i.e. -p[i][k]) */</font>
572 00564 <font class="keywordflow">for</font> (k=0; k&lt;i; k++)
573 00565 <font class="keywordflow">if</font> (p[k][i] == p[i][j]) p[k][j] = p[k][i];
574 00566
575 00567 <font class="comment">/* then, i&lt;k&lt;j so look for p[i][j] == -p[i][k] */</font>
576 00568 <font class="keywordflow">for</font> (k=i+1; k&lt;j; k++)
577 00569 <font class="keywordflow">if</font> (p[i][k] == -p[i][j]) p[k][j] = -p[i][k];
578 00570
579 00571 <font class="comment">/* last, k&gt;j same search but */</font>
580 00572 <font class="keywordflow">for</font> (k=j+1; k&lt;n; k++)
581 00573 <font class="keywordflow">if</font> (p[i][k] == -p[i][j]) p[j][k] = p[i][k];
582 00574 }
583 00575 }
584 00576 }
585 00577
586 00578 <font class="comment">/* walk thru the dag to compute the partial order (and optionally</font>
587 00579 <font class="comment"> the permutation)</font>
588 00580 <font class="comment"> Note: this is not the fastest way to do it but it takes</font>
589 00581 <font class="comment"> negligible time compared to a single call of PolyhedronLTQ !</font>
590 00582 <font class="comment"> Each macro-step (while loop) gives the same logical time to all</font>
591 00583 <font class="comment"> found roots and optionally add these nodes in the permutation vector</font>
592 00584 <font class="comment"> */</font>
593 00585
594 00586 t = 0; <font class="comment">/* current logical time, assigned to current roots and</font>
595 00587 <font class="comment"> increased by 1 at the end of each step */</font>
596 00588 nb = 0; <font class="comment">/* number of processed nodes (have been given a time) */</font>
597 00589 <font class="keywordflow">while</font> (nb&lt;n) {
598 00590 <font class="keywordflow">for</font> (i=0; i&lt;n; i++) { <font class="comment">/* search for roots */</font>
599 00591 <font class="comment">/* for any node, if it's not already been given a logical time</font>
600 00592 <font class="comment"> then walk thru the node row; if we find a 1 at some column j,</font>
601 00593 <font class="comment"> it means that node j preceeds the current node, so it is not</font>
602 00594 <font class="comment"> a root */</font>
603 00595 <font class="keywordflow">if</font> (time[i]&lt;0) {
604 00596 isroot = 1; <font class="comment">/* assume that it is until it is definitely not */</font>
605 00597 <font class="comment">/* first search on a column */</font>
606 00598 <font class="keywordflow">for</font> (j=0; j&lt;i; j++) {
607 00599 <font class="keywordflow">if</font> (p[j][i]==-1) { <font class="comment">/* found a node lower than it */</font>
608 00600 isroot = 0; <font class="keywordflow">break</font>;
609 00601 }
610 00602 }
611 00603 <font class="keywordflow">if</font> <font class="comment">/*still*/</font> (isroot)
612 00604 <font class="keywordflow">for</font> (j=i+1; j&lt;n; j++) {
613 00605 <font class="keywordflow">if</font> (p[i][j]==1) { <font class="comment">/* greater than this node */</font>
614 00606 isroot = 0; <font class="keywordflow">break</font>;
615 00607 }
616 00608 }
617 00609 <font class="keywordflow">if</font> (isroot) { <font class="comment">/* then it definitely is */</font>
618 00610 time[i] = t; <font class="comment">/* this node gets current logical time */</font>
619 00611 <font class="keywordflow">if</font> (pvect)
620 00612 pvect[nb] = i+1; <font class="comment">/* nodes will be numbered from 1 to n */</font>
621 00613 nb++; <font class="comment">/* one more node processed */</font>
622 00614 }
623 00615 }
624 00616 }
625 00617 <font class="comment">/* now make roots become neutral, i.e. non comparable with all other nodes */</font>
626 00618 <font class="keywordflow">for</font> (i=0; i&lt;n; i++) {
627 00619 <font class="keywordflow">if</font> (time[i]==t) {
628 00620 <font class="keywordflow">for</font> (j=0; j&lt;i; j++)
629 00621 p[j][i] = 0;
630 00622 <font class="keywordflow">for</font> (j=i+1; j&lt;n; j++)
631 00623 p[i][j] = 0;
632 00624 }
633 00625 }
634 00626 t++; <font class="comment">/* ready for next set of root nodes */</font>
635 00627 }
636 00628
637 00629 free (p); <font class="comment">/* let's be clean, collect the garbage */</font>
638 00630 free (dag);
639 00631 <font class="keywordflow">return</font> 1;
640 00632 } <font class="comment">/* PolyhedronTSort */</font>
641 00633
642 00634
643 00635
644 00636
645 </pre></div><hr><address align="right"><small>Generated on Fri Nov 8 12:10:06 2002 for Polylib by
646 <a href="http://www.doxygen.org/index.html">
647 <img src="doxygen.png" alt="doxygen" align="middle" border=0
648 width=110 height=53></a>1.2.15 </small></address>
649 </body>
650 </html>
Polyhedral library