1 /****************************************************************************************/
2 /* */
3 /* vinci_computation.c */
4 /* */
5 /****************************************************************************************/
6 /* */
7 /* Authors: Benno Bueeler (bueeler@ifor.math.ethz.ch) */
8 /* and */
9 /* Andreas Enge (enge@ifor.math.ethz.ch) */
10 /* Institute for Operations Research */
11 /* Swiss Federal Institute of Technology Zurich */
12 /* Switzerland */
13 /* */
14 /* Last Changes: February 2, 2001 */
15 /* */
16 /****************************************************************************************/
17 /* */
18 /* the main computational routines of the package */
19 /* */
20 /****************************************************************************************/
21
22 #include "vinci.h"
23
24 /****************************************************************************************/
25
factorial(int n)26 rational factorial (int n)
27 /* calculates n! by calculating all i! for 0 <= i <= G_d and storing them in an */
28 /* array */
29
30 { static rational *fact;
31 static boolean first_call = TRUE;
32 int i;
33
34 if (first_call)
35 {
36 fact = create_fact ();
37 fact [0] = 1;
38 for (i = 1; i <= G_d; i++) fact [i] = i * fact [i - 1];
39 first_call = FALSE;
40 }
41
42 return (fact [n]);
43 }
44
45 /****************************************************************************************/
46 /* computations with simplices and determinants */
47 /****************************************************************************************/
48
det_and_invert(rational ** A,int rows,int columns,boolean verbose)49 rational det_and_invert (rational **A, int rows, int columns, boolean verbose)
50 /* Let B be the quadratic submatrix of A defined by the first "rows" rows and */
51 /* columns of A. Then this function turns B into an upper triangular matrix by means */
52 /* of Gaussian pivot operations. At the same time, the colums "rows + 1" up to */
53 /* "columns" of A are replaced by B^-1 times these columns. As a byproduct, the de- */
54 /* terminant of B (modulo its sign) is returned. */
55 /* If verbose is TRUE an error message is output for zero volume. */
56
57 { rational dummy, det = 1, *row_pointer;
58 int i, j, k, r;
59 #if PIVOTING == 2
60 int s;
61 #elif PIVOTING == 0
62 boolean pivot_found = FALSE;
63 #endif
64
65 for (k = 0; k < rows; k++)
66 { /* look for a pivot element */
67 r = k;
68
69 #if PIVOTING == 0
70 /* choose the first row r with r >= k and A (r, k) > MIN_PIVOT */
71 for (i = k + 1; i < rows && !pivot_found; i++)
72 { dummy = fabs (A [i] [k]);
73 if (dummy > MIN_PIVOT)
74 { r = i;
75 pivot_found = TRUE;
76 }
77 else if (dummy > fabs (A [r] [k])) r = i;
78 }
79 #elif PIVOTING == 2
80 /* total pivoting; choose row r and column s with maximal |A (r, s)|, but only if */
81 /* columns is equal to rows (no equations have to be solved); otherwise use par- */
82 /* tial pivoting */
83 if (columns == rows)
84 for (i = k; i < G_d; i++)
85 for (j = k; j < G_d; j++)
86 if (fabs (A [i] [j]) > fabs (A [r] [s]))
87 { r = i;
88 s = j;}
89 else
90 for (i = k + 1; i < rows; i++)
91 if (fabs (A [i] [k]) > fabs (A [r] [k])) r = i;
92 #else
93 /* partial pivoting; choose row r with maximal |A (r, k)| */
94 for (i = k + 1; i < rows; i++)
95 if (fabs (A [i] [k]) > fabs (A [r] [k])) r = i;
96 #endif
97
98 if (r != k) /* exchange the rows r and k */
99 { row_pointer = A [k];
100 A [k] = A [r];
101 A [r] = row_pointer;
102 }
103
104 /* check for singularity */
105 if (fabs (A [k] [k]) < EPSILON)
106 { if (verbose)
107 { fprintf (stderr, "\n***** ERROR: (Almost) singular matrix in 'det_and_");
108 fprintf (stderr, "invert';\n One pivot is %20.18e.\n", (real) A [k] [k]);
109 }
110 return 0;
111 }
112
113 /* now make an elimination step */
114 for (i = k + 1; i < rows; i++)
115 { dummy = A [i] [k] / A [k] [k]; /* store the factor once for the row */
116 for (j = k + 1; j < columns; j++) A [i] [j] -= dummy * A [k] [j];
117 };
118 det *= A [k] [k];
119 }
120
121 /* now B is an upper triangular matrix; calculate B^-1 times the last columns */
122 for (k = rows; k < columns; k++)
123 for (i = rows - 1; i >= 0; i--)
124 { /* compute (B^-1 * A_k)_i */
125 for (j = i + 1; j < rows; j++) A [i] [k] -= A [i] [j] * A [j] [k];
126 A [i] [k] /= A [i] [i];
127 }
128
129 return det;
130 }
131
132 /****************************************************************************************/
133
simplex_volume(T_VertexSet S,rational * volume,boolean verbose)134 void simplex_volume (T_VertexSet S, rational *volume, boolean verbose)
135 /* From a set of vertices which are presumed of number G_d + 1 and affinely indepen- */
136 /* dent, the function computes G_d! times the volume of the corresponding simplex */
137 /* and stores it in volume. */
138 /* If verbose is TRUE an error message is output for zero volume. */
139
140 { static rational **A;
141 static boolean first_call = TRUE;
142 rational dummy;
143 int i, j;
144
145 if (first_call)
146 {
147 A = create_matrix (G_d, G_d);
148 first_call = FALSE;
149 }
150
151 /* copy the relevant information into A */
152 for (j = 0; j < G_d; j++)
153 { dummy = S.loe [G_d] -> coords [j];
154 for (i = 0; i < G_d; i++)
155 A [i] [j] = (S.loe [i] -> coords [j]) - dummy;
156 }
157
158 *volume = det_and_invert (A, G_d, G_d, verbose);
159
160 if (*volume < 0) *volume = - (*volume);
161
162 #ifdef STATISTICS
163 if (*volume > EPSILON) update_statistics ((real) *volume / factorial (G_d));
164 #endif
165
166 }
167
168 /****************************************************************************************/
169 /* routines for computing orthonormal bases */
170 /****************************************************************************************/
171
add_orthonormal(int d,T_VertexSet face,rational ** H,T_Vertex * vertex)172 rational add_orthonormal (int d, T_VertexSet face, rational **H, T_Vertex *vertex)
173 /* The function assumes Householder vectors of the linear subspace associated with */
174 /* the variable face in the first d-1 rows of H. It computes the Householder vector */
175 /* for the new vertex minus a vertex of face and stores it into the dth row of H. */
176 /* The first k-1 zero entries of the Householder vector k are neglected, its essen- */
177 /* tial part is stored in the columns k to dimension. All Householder vectors are */
178 /* normalised. This may be changed, but then do not forget to divide by its norm in */
179 /* later steps! */
180 /* Return value is the distance of the vertex to face. */
181
182 { int i, j;
183 rational scalar_product, alpha_squared, alpha, divisor, distance, mindistance = 1e99;
184
185 /* look for the vertex in face with the smallest distance to the given vertex */
186 for (i = 0; i <= face.lastel; i++)
187 {
188 distance = 0;
189 for (j = 0; j < G_d; j++)
190 distance += (vertex -> coords [j] - face.loe [i] -> coords [j])
191 * (vertex -> coords [j] - face.loe [i] -> coords [j]);
192 if (distance < mindistance)
193 { mindistance = distance;
194 for (j = 0; j < G_d; j++)
195 H [d-1] [j] = vertex -> coords [j] - face.loe [i] -> coords [j];
196 }
197 }
198
199 /* multiply from the left with the previous Householder matrices */
200 for (i = 0; i < d-1; i++)
201 { scalar_product = 0;
202 for (j = i; j < G_d; j++)
203 scalar_product += H [i] [j] * H [d-1] [j];
204 scalar_product *= 2;
205 for (j = i; j < G_d; j++)
206 H [d-1] [j] -= scalar_product * H [i] [j];
207 }
208
209 /* compute the new Householder vector and the distance */
210 alpha_squared = 0;
211 for (j = d-1; j < G_d; j++)
212 alpha_squared += H [d-1] [j] * H [d-1] [j];
213 alpha = sqrt (alpha_squared);
214 if (H [d-1] [d-1] < 0)
215 alpha = - alpha;
216
217 divisor = sqrt (2 * (alpha_squared + alpha * H [d-1] [d-1]));
218 H [d-1] [d-1] += alpha;
219 for (j = d-1; j < G_d; j++)
220 H [d-1] [j] /= divisor;
221
222 return fabs (alpha);
223
224 }
225
226 /****************************************************************************************/
227
orthonormal(int d,T_VertexSet face,rational ** H)228 rational orthonormal (int d, T_VertexSet face, rational **H)
229 /* The function computes Householder vectors for the linear subspace associated with */
230 /* the variable face. The vectors are returned via the first d rows of H. */
231 /* This routine uses pivoting techniques and should be numerically stable. */
232 /* The return value is d! times the volume of the orthonormalised simplex. */
233
234 { static rational **local_H;
235 static int m = 0;
236 /* the number of rows in local_H */
237 rational *dummy_row;
238 int i, j, k, maxindex = 0;
239 rational scalar_product, alpha_squared = -1, alpha, divisor;
240 rational volume = 1;
241
242 /* create local_H in the correct dimension */
243 if (m == 0)
244 {
245 m = face.lastel;
246 local_H = create_matrix (m, G_d + 1);
247 /* The last component of each row will contain the norm of the essential part */
248 /* of the vector */
249 }
250 else if (face.lastel > m)
251 {
252 redim_matrix (&local_H, m, face.lastel, G_d + 1);
253 m = face.lastel;
254 }
255
256 /* copy the spanning vectors into local_H; the zeroth vertex of face [d] is trans- */
257 /* lated into the origin */
258 for (i = 0; i < face.lastel; i++)
259 { local_H [i] [G_d] = 0;
260 for (j = 0; j < G_d; j++)
261 { local_H [i] [j] = face.loe [i+1] -> coords [j]
262 - face.loe [0 ] -> coords [j];
263 local_H [i] [G_d] += local_H [i] [j] * local_H [i] [j];
264 }
265 if (local_H [i] [G_d] > alpha_squared)
266 { alpha_squared = local_H [i] [G_d];
267 maxindex = i;
268 }
269 }
270
271 for (k = 0; k < d; k++)
272 { /* compute k-th Householder vector; choose the element whose relevant part has */
273 /* the biggest norm as next element */
274
275 /* If alpha_squared is 0 matrix has not got the full rank. Stop orthonormalising. */
276 if (alpha_squared / volume < EPSILON)
277 { return 0;
278 }
279
280 /* otherwise exchange rows */
281 dummy_row = local_H [k];
282 local_H [k] = local_H [maxindex];
283 local_H [maxindex] = dummy_row;
284
285 /* compute the new Householder vector */
286 alpha = sqrt (alpha_squared);
287 volume *= alpha;
288 if (local_H [k] [k] < 0)
289 alpha = - alpha;
290 divisor = sqrt (2 * (alpha_squared + alpha * local_H [k] [k]));
291 local_H [k] [k] += alpha;
292 for (j = k; j < G_d; j++)
293 local_H [k] [j] /= divisor;
294
295 /* apply the Householder matrix to the resting rows of local_H; at the same time */
296 /* compute the norms of the relevant parts */
297 alpha_squared = -1;
298 for (i = k + 1; i < face.lastel; i++)
299 {
300 scalar_product = 0;
301 for (j = k; j < G_d; j++)
302 scalar_product += local_H [k] [j] * local_H [i] [j];
303 scalar_product *= 2;
304 for (j = k; j < G_d; j++)
305 local_H [i] [j] -= scalar_product * local_H [k] [j];
306
307 local_H [i] [G_d] -= local_H [i] [k] * local_H [i] [k];
308 if (local_H [i] [G_d] > alpha_squared)
309 { alpha_squared = local_H [i] [G_d];
310 maxindex = i;
311 }
312 }
313
314 } /* for */
315
316 /* copy the results into H */
317 for (i = 0; i < d; i++)
318 for (j = 0; j < G_d; j++)
319 H [i] [j] = local_H [i] [j];
320
321 return volume;
322 }
323
324 /****************************************************************************************/
325