|
0
|
1 /* Integer matrix math routines
|
|
|
2 Copyright (C) 2003, 2004, 2005, 2007, 2008 Free Software Foundation, Inc.
|
|
|
3 Contributed by Daniel Berlin <dberlin@dberlin.org>.
|
|
|
4
|
|
|
5 This file is part of GCC.
|
|
|
6
|
|
|
7 GCC is free software; you can redistribute it and/or modify it under
|
|
|
8 the terms of the GNU General Public License as published by the Free
|
|
|
9 Software Foundation; either version 3, or (at your option) any later
|
|
|
10 version.
|
|
|
11
|
|
|
12 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
|
13 WARRANTY; without even the implied warranty of MERCHANTABILITY or
|
|
|
14 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
|
|
|
15 for more details.
|
|
|
16
|
|
|
17 You should have received a copy of the GNU General Public License
|
|
|
18 along with GCC; see the file COPYING3. If not see
|
|
|
19 <http://www.gnu.org/licenses/>. */
|
|
|
20
|
|
|
21 #include "config.h"
|
|
|
22 #include "system.h"
|
|
|
23 #include "coretypes.h"
|
|
|
24 #include "tm.h"
|
|
|
25 #include "ggc.h"
|
|
|
26 #include "tree.h"
|
|
|
27 #include "tree-flow.h"
|
|
|
28 #include "lambda.h"
|
|
|
29
|
|
|
30 static void lambda_matrix_get_column (lambda_matrix, int, int,
|
|
|
31 lambda_vector);
|
|
|
32
|
|
|
33 /* Allocate a matrix of M rows x N cols. */
|
|
|
34
|
|
|
35 lambda_matrix
|
|
|
36 lambda_matrix_new (int m, int n)
|
|
|
37 {
|
|
|
38 lambda_matrix mat;
|
|
|
39 int i;
|
|
|
40
|
|
|
41 mat = GGC_NEWVEC (lambda_vector, m);
|
|
|
42
|
|
|
43 for (i = 0; i < m; i++)
|
|
|
44 mat[i] = lambda_vector_new (n);
|
|
|
45
|
|
|
46 return mat;
|
|
|
47 }
|
|
|
48
|
|
|
49 /* Copy the elements of M x N matrix MAT1 to MAT2. */
|
|
|
50
|
|
|
51 void
|
|
|
52 lambda_matrix_copy (lambda_matrix mat1, lambda_matrix mat2,
|
|
|
53 int m, int n)
|
|
|
54 {
|
|
|
55 int i;
|
|
|
56
|
|
|
57 for (i = 0; i < m; i++)
|
|
|
58 lambda_vector_copy (mat1[i], mat2[i], n);
|
|
|
59 }
|
|
|
60
|
|
|
61 /* Store the N x N identity matrix in MAT. */
|
|
|
62
|
|
|
63 void
|
|
|
64 lambda_matrix_id (lambda_matrix mat, int size)
|
|
|
65 {
|
|
|
66 int i, j;
|
|
|
67
|
|
|
68 for (i = 0; i < size; i++)
|
|
|
69 for (j = 0; j < size; j++)
|
|
|
70 mat[i][j] = (i == j) ? 1 : 0;
|
|
|
71 }
|
|
|
72
|
|
|
73 /* Return true if MAT is the identity matrix of SIZE */
|
|
|
74
|
|
|
75 bool
|
|
|
76 lambda_matrix_id_p (lambda_matrix mat, int size)
|
|
|
77 {
|
|
|
78 int i, j;
|
|
|
79 for (i = 0; i < size; i++)
|
|
|
80 for (j = 0; j < size; j++)
|
|
|
81 {
|
|
|
82 if (i == j)
|
|
|
83 {
|
|
|
84 if (mat[i][j] != 1)
|
|
|
85 return false;
|
|
|
86 }
|
|
|
87 else
|
|
|
88 {
|
|
|
89 if (mat[i][j] != 0)
|
|
|
90 return false;
|
|
|
91 }
|
|
|
92 }
|
|
|
93 return true;
|
|
|
94 }
|
|
|
95
|
|
|
96 /* Negate the elements of the M x N matrix MAT1 and store it in MAT2. */
|
|
|
97
|
|
|
98 void
|
|
|
99 lambda_matrix_negate (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
|
|
|
100 {
|
|
|
101 int i;
|
|
|
102
|
|
|
103 for (i = 0; i < m; i++)
|
|
|
104 lambda_vector_negate (mat1[i], mat2[i], n);
|
|
|
105 }
|
|
|
106
|
|
|
107 /* Take the transpose of matrix MAT1 and store it in MAT2.
|
|
|
108 MAT1 is an M x N matrix, so MAT2 must be N x M. */
|
|
|
109
|
|
|
110 void
|
|
|
111 lambda_matrix_transpose (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
|
|
|
112 {
|
|
|
113 int i, j;
|
|
|
114
|
|
|
115 for (i = 0; i < n; i++)
|
|
|
116 for (j = 0; j < m; j++)
|
|
|
117 mat2[i][j] = mat1[j][i];
|
|
|
118 }
|
|
|
119
|
|
|
120
|
|
|
121 /* Add two M x N matrices together: MAT3 = MAT1+MAT2. */
|
|
|
122
|
|
|
123 void
|
|
|
124 lambda_matrix_add (lambda_matrix mat1, lambda_matrix mat2,
|
|
|
125 lambda_matrix mat3, int m, int n)
|
|
|
126 {
|
|
|
127 int i;
|
|
|
128
|
|
|
129 for (i = 0; i < m; i++)
|
|
|
130 lambda_vector_add (mat1[i], mat2[i], mat3[i], n);
|
|
|
131 }
|
|
|
132
|
|
|
133 /* MAT3 = CONST1 * MAT1 + CONST2 * MAT2. All matrices are M x N. */
|
|
|
134
|
|
|
135 void
|
|
|
136 lambda_matrix_add_mc (lambda_matrix mat1, int const1,
|
|
|
137 lambda_matrix mat2, int const2,
|
|
|
138 lambda_matrix mat3, int m, int n)
|
|
|
139 {
|
|
|
140 int i;
|
|
|
141
|
|
|
142 for (i = 0; i < m; i++)
|
|
|
143 lambda_vector_add_mc (mat1[i], const1, mat2[i], const2, mat3[i], n);
|
|
|
144 }
|
|
|
145
|
|
|
146 /* Multiply two matrices: MAT3 = MAT1 * MAT2.
|
|
|
147 MAT1 is an M x R matrix, and MAT2 is R x N. The resulting MAT2
|
|
|
148 must therefore be M x N. */
|
|
|
149
|
|
|
150 void
|
|
|
151 lambda_matrix_mult (lambda_matrix mat1, lambda_matrix mat2,
|
|
|
152 lambda_matrix mat3, int m, int r, int n)
|
|
|
153 {
|
|
|
154
|
|
|
155 int i, j, k;
|
|
|
156
|
|
|
157 for (i = 0; i < m; i++)
|
|
|
158 {
|
|
|
159 for (j = 0; j < n; j++)
|
|
|
160 {
|
|
|
161 mat3[i][j] = 0;
|
|
|
162 for (k = 0; k < r; k++)
|
|
|
163 mat3[i][j] += mat1[i][k] * mat2[k][j];
|
|
|
164 }
|
|
|
165 }
|
|
|
166 }
|
|
|
167
|
|
|
168 /* Get column COL from the matrix MAT and store it in VEC. MAT has
|
|
|
169 N rows, so the length of VEC must be N. */
|
|
|
170
|
|
|
171 static void
|
|
|
172 lambda_matrix_get_column (lambda_matrix mat, int n, int col,
|
|
|
173 lambda_vector vec)
|
|
|
174 {
|
|
|
175 int i;
|
|
|
176
|
|
|
177 for (i = 0; i < n; i++)
|
|
|
178 vec[i] = mat[i][col];
|
|
|
179 }
|
|
|
180
|
|
|
181 /* Delete rows r1 to r2 (not including r2). */
|
|
|
182
|
|
|
183 void
|
|
|
184 lambda_matrix_delete_rows (lambda_matrix mat, int rows, int from, int to)
|
|
|
185 {
|
|
|
186 int i;
|
|
|
187 int dist;
|
|
|
188 dist = to - from;
|
|
|
189
|
|
|
190 for (i = to; i < rows; i++)
|
|
|
191 mat[i - dist] = mat[i];
|
|
|
192
|
|
|
193 for (i = rows - dist; i < rows; i++)
|
|
|
194 mat[i] = NULL;
|
|
|
195 }
|
|
|
196
|
|
|
197 /* Swap rows R1 and R2 in matrix MAT. */
|
|
|
198
|
|
|
199 void
|
|
|
200 lambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2)
|
|
|
201 {
|
|
|
202 lambda_vector row;
|
|
|
203
|
|
|
204 row = mat[r1];
|
|
|
205 mat[r1] = mat[r2];
|
|
|
206 mat[r2] = row;
|
|
|
207 }
|
|
|
208
|
|
|
209 /* Add a multiple of row R1 of matrix MAT with N columns to row R2:
|
|
|
210 R2 = R2 + CONST1 * R1. */
|
|
|
211
|
|
|
212 void
|
|
|
213 lambda_matrix_row_add (lambda_matrix mat, int n, int r1, int r2, int const1)
|
|
|
214 {
|
|
|
215 int i;
|
|
|
216
|
|
|
217 if (const1 == 0)
|
|
|
218 return;
|
|
|
219
|
|
|
220 for (i = 0; i < n; i++)
|
|
|
221 mat[r2][i] += const1 * mat[r1][i];
|
|
|
222 }
|
|
|
223
|
|
|
224 /* Negate row R1 of matrix MAT which has N columns. */
|
|
|
225
|
|
|
226 void
|
|
|
227 lambda_matrix_row_negate (lambda_matrix mat, int n, int r1)
|
|
|
228 {
|
|
|
229 lambda_vector_negate (mat[r1], mat[r1], n);
|
|
|
230 }
|
|
|
231
|
|
|
232 /* Multiply row R1 of matrix MAT with N columns by CONST1. */
|
|
|
233
|
|
|
234 void
|
|
|
235 lambda_matrix_row_mc (lambda_matrix mat, int n, int r1, int const1)
|
|
|
236 {
|
|
|
237 int i;
|
|
|
238
|
|
|
239 for (i = 0; i < n; i++)
|
|
|
240 mat[r1][i] *= const1;
|
|
|
241 }
|
|
|
242
|
|
|
243 /* Exchange COL1 and COL2 in matrix MAT. M is the number of rows. */
|
|
|
244
|
|
|
245 void
|
|
|
246 lambda_matrix_col_exchange (lambda_matrix mat, int m, int col1, int col2)
|
|
|
247 {
|
|
|
248 int i;
|
|
|
249 int tmp;
|
|
|
250 for (i = 0; i < m; i++)
|
|
|
251 {
|
|
|
252 tmp = mat[i][col1];
|
|
|
253 mat[i][col1] = mat[i][col2];
|
|
|
254 mat[i][col2] = tmp;
|
|
|
255 }
|
|
|
256 }
|
|
|
257
|
|
|
258 /* Add a multiple of column C1 of matrix MAT with M rows to column C2:
|
|
|
259 C2 = C2 + CONST1 * C1. */
|
|
|
260
|
|
|
261 void
|
|
|
262 lambda_matrix_col_add (lambda_matrix mat, int m, int c1, int c2, int const1)
|
|
|
263 {
|
|
|
264 int i;
|
|
|
265
|
|
|
266 if (const1 == 0)
|
|
|
267 return;
|
|
|
268
|
|
|
269 for (i = 0; i < m; i++)
|
|
|
270 mat[i][c2] += const1 * mat[i][c1];
|
|
|
271 }
|
|
|
272
|
|
|
273 /* Negate column C1 of matrix MAT which has M rows. */
|
|
|
274
|
|
|
275 void
|
|
|
276 lambda_matrix_col_negate (lambda_matrix mat, int m, int c1)
|
|
|
277 {
|
|
|
278 int i;
|
|
|
279
|
|
|
280 for (i = 0; i < m; i++)
|
|
|
281 mat[i][c1] *= -1;
|
|
|
282 }
|
|
|
283
|
|
|
284 /* Multiply column C1 of matrix MAT with M rows by CONST1. */
|
|
|
285
|
|
|
286 void
|
|
|
287 lambda_matrix_col_mc (lambda_matrix mat, int m, int c1, int const1)
|
|
|
288 {
|
|
|
289 int i;
|
|
|
290
|
|
|
291 for (i = 0; i < m; i++)
|
|
|
292 mat[i][c1] *= const1;
|
|
|
293 }
|
|
|
294
|
|
|
295 /* Compute the inverse of the N x N matrix MAT and store it in INV.
|
|
|
296
|
|
|
297 We don't _really_ compute the inverse of MAT. Instead we compute
|
|
|
298 det(MAT)*inv(MAT), and we return det(MAT) to the caller as the function
|
|
|
299 result. This is necessary to preserve accuracy, because we are dealing
|
|
|
300 with integer matrices here.
|
|
|
301
|
|
|
302 The algorithm used here is a column based Gauss-Jordan elimination on MAT
|
|
|
303 and the identity matrix in parallel. The inverse is the result of applying
|
|
|
304 the same operations on the identity matrix that reduce MAT to the identity
|
|
|
305 matrix.
|
|
|
306
|
|
|
307 When MAT is a 2 x 2 matrix, we don't go through the whole process, because
|
|
|
308 it is easily inverted by inspection and it is a very common case. */
|
|
|
309
|
|
|
310 static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int);
|
|
|
311
|
|
|
312 int
|
|
|
313 lambda_matrix_inverse (lambda_matrix mat, lambda_matrix inv, int n)
|
|
|
314 {
|
|
|
315 if (n == 2)
|
|
|
316 {
|
|
|
317 int a, b, c, d, det;
|
|
|
318 a = mat[0][0];
|
|
|
319 b = mat[1][0];
|
|
|
320 c = mat[0][1];
|
|
|
321 d = mat[1][1];
|
|
|
322 inv[0][0] = d;
|
|
|
323 inv[0][1] = -c;
|
|
|
324 inv[1][0] = -b;
|
|
|
325 inv[1][1] = a;
|
|
|
326 det = (a * d - b * c);
|
|
|
327 if (det < 0)
|
|
|
328 {
|
|
|
329 det *= -1;
|
|
|
330 inv[0][0] *= -1;
|
|
|
331 inv[1][0] *= -1;
|
|
|
332 inv[0][1] *= -1;
|
|
|
333 inv[1][1] *= -1;
|
|
|
334 }
|
|
|
335 return det;
|
|
|
336 }
|
|
|
337 else
|
|
|
338 return lambda_matrix_inverse_hard (mat, inv, n);
|
|
|
339 }
|
|
|
340
|
|
|
341 /* If MAT is not a special case, invert it the hard way. */
|
|
|
342
|
|
|
343 static int
|
|
|
344 lambda_matrix_inverse_hard (lambda_matrix mat, lambda_matrix inv, int n)
|
|
|
345 {
|
|
|
346 lambda_vector row;
|
|
|
347 lambda_matrix temp;
|
|
|
348 int i, j;
|
|
|
349 int determinant;
|
|
|
350
|
|
|
351 temp = lambda_matrix_new (n, n);
|
|
|
352 lambda_matrix_copy (mat, temp, n, n);
|
|
|
353 lambda_matrix_id (inv, n);
|
|
|
354
|
|
|
355 /* Reduce TEMP to a lower triangular form, applying the same operations on
|
|
|
356 INV which starts as the identity matrix. N is the number of rows and
|
|
|
357 columns. */
|
|
|
358 for (j = 0; j < n; j++)
|
|
|
359 {
|
|
|
360 row = temp[j];
|
|
|
361
|
|
|
362 /* Make every element in the current row positive. */
|
|
|
363 for (i = j; i < n; i++)
|
|
|
364 if (row[i] < 0)
|
|
|
365 {
|
|
|
366 lambda_matrix_col_negate (temp, n, i);
|
|
|
367 lambda_matrix_col_negate (inv, n, i);
|
|
|
368 }
|
|
|
369
|
|
|
370 /* Sweep the upper triangle. Stop when only the diagonal element in the
|
|
|
371 current row is nonzero. */
|
|
|
372 while (lambda_vector_first_nz (row, n, j + 1) < n)
|
|
|
373 {
|
|
|
374 int min_col = lambda_vector_min_nz (row, n, j);
|
|
|
375 lambda_matrix_col_exchange (temp, n, j, min_col);
|
|
|
376 lambda_matrix_col_exchange (inv, n, j, min_col);
|
|
|
377
|
|
|
378 for (i = j + 1; i < n; i++)
|
|
|
379 {
|
|
|
380 int factor;
|
|
|
381
|
|
|
382 factor = -1 * row[i];
|
|
|
383 if (row[j] != 1)
|
|
|
384 factor /= row[j];
|
|
|
385
|
|
|
386 lambda_matrix_col_add (temp, n, j, i, factor);
|
|
|
387 lambda_matrix_col_add (inv, n, j, i, factor);
|
|
|
388 }
|
|
|
389 }
|
|
|
390 }
|
|
|
391
|
|
|
392 /* Reduce TEMP from a lower triangular to the identity matrix. Also compute
|
|
|
393 the determinant, which now is simply the product of the elements on the
|
|
|
394 diagonal of TEMP. If one of these elements is 0, the matrix has 0 as an
|
|
|
395 eigenvalue so it is singular and hence not invertible. */
|
|
|
396 determinant = 1;
|
|
|
397 for (j = n - 1; j >= 0; j--)
|
|
|
398 {
|
|
|
399 int diagonal;
|
|
|
400
|
|
|
401 row = temp[j];
|
|
|
402 diagonal = row[j];
|
|
|
403
|
|
|
404 /* The matrix must not be singular. */
|
|
|
405 gcc_assert (diagonal);
|
|
|
406
|
|
|
407 determinant = determinant * diagonal;
|
|
|
408
|
|
|
409 /* If the diagonal is not 1, then multiply the each row by the
|
|
|
410 diagonal so that the middle number is now 1, rather than a
|
|
|
411 rational. */
|
|
|
412 if (diagonal != 1)
|
|
|
413 {
|
|
|
414 for (i = 0; i < j; i++)
|
|
|
415 lambda_matrix_col_mc (inv, n, i, diagonal);
|
|
|
416 for (i = j + 1; i < n; i++)
|
|
|
417 lambda_matrix_col_mc (inv, n, i, diagonal);
|
|
|
418
|
|
|
419 row[j] = diagonal = 1;
|
|
|
420 }
|
|
|
421
|
|
|
422 /* Sweep the lower triangle column wise. */
|
|
|
423 for (i = j - 1; i >= 0; i--)
|
|
|
424 {
|
|
|
425 if (row[i])
|
|
|
426 {
|
|
|
427 int factor = -row[i];
|
|
|
428 lambda_matrix_col_add (temp, n, j, i, factor);
|
|
|
429 lambda_matrix_col_add (inv, n, j, i, factor);
|
|
|
430 }
|
|
|
431
|
|
|
432 }
|
|
|
433 }
|
|
|
434
|
|
|
435 return determinant;
|
|
|
436 }
|
|
|
437
|
|
|
438 /* Decompose a N x N matrix MAT to a product of a lower triangular H
|
|
|
439 and a unimodular U matrix such that MAT = H.U. N is the size of
|
|
|
440 the rows of MAT. */
|
|
|
441
|
|
|
442 void
|
|
|
443 lambda_matrix_hermite (lambda_matrix mat, int n,
|
|
|
444 lambda_matrix H, lambda_matrix U)
|
|
|
445 {
|
|
|
446 lambda_vector row;
|
|
|
447 int i, j, factor, minimum_col;
|
|
|
448
|
|
|
449 lambda_matrix_copy (mat, H, n, n);
|
|
|
450 lambda_matrix_id (U, n);
|
|
|
451
|
|
|
452 for (j = 0; j < n; j++)
|
|
|
453 {
|
|
|
454 row = H[j];
|
|
|
455
|
|
|
456 /* Make every element of H[j][j..n] positive. */
|
|
|
457 for (i = j; i < n; i++)
|
|
|
458 {
|
|
|
459 if (row[i] < 0)
|
|
|
460 {
|
|
|
461 lambda_matrix_col_negate (H, n, i);
|
|
|
462 lambda_vector_negate (U[i], U[i], n);
|
|
|
463 }
|
|
|
464 }
|
|
|
465
|
|
|
466 /* Stop when only the diagonal element is nonzero. */
|
|
|
467 while (lambda_vector_first_nz (row, n, j + 1) < n)
|
|
|
468 {
|
|
|
469 minimum_col = lambda_vector_min_nz (row, n, j);
|
|
|
470 lambda_matrix_col_exchange (H, n, j, minimum_col);
|
|
|
471 lambda_matrix_row_exchange (U, j, minimum_col);
|
|
|
472
|
|
|
473 for (i = j + 1; i < n; i++)
|
|
|
474 {
|
|
|
475 factor = row[i] / row[j];
|
|
|
476 lambda_matrix_col_add (H, n, j, i, -1 * factor);
|
|
|
477 lambda_matrix_row_add (U, n, i, j, factor);
|
|
|
478 }
|
|
|
479 }
|
|
|
480 }
|
|
|
481 }
|
|
|
482
|
|
|
483 /* Given an M x N integer matrix A, this function determines an M x
|
|
|
484 M unimodular matrix U, and an M x N echelon matrix S such that
|
|
|
485 "U.A = S". This decomposition is also known as "right Hermite".
|
|
|
486
|
|
|
487 Ref: Algorithm 2.1 page 33 in "Loop Transformations for
|
|
|
488 Restructuring Compilers" Utpal Banerjee. */
|
|
|
489
|
|
|
490 void
|
|
|
491 lambda_matrix_right_hermite (lambda_matrix A, int m, int n,
|
|
|
492 lambda_matrix S, lambda_matrix U)
|
|
|
493 {
|
|
|
494 int i, j, i0 = 0;
|
|
|
495
|
|
|
496 lambda_matrix_copy (A, S, m, n);
|
|
|
497 lambda_matrix_id (U, m);
|
|
|
498
|
|
|
499 for (j = 0; j < n; j++)
|
|
|
500 {
|
|
|
501 if (lambda_vector_first_nz (S[j], m, i0) < m)
|
|
|
502 {
|
|
|
503 ++i0;
|
|
|
504 for (i = m - 1; i >= i0; i--)
|
|
|
505 {
|
|
|
506 while (S[i][j] != 0)
|
|
|
507 {
|
|
|
508 int sigma, factor, a, b;
|
|
|
509
|
|
|
510 a = S[i-1][j];
|
|
|
511 b = S[i][j];
|
|
|
512 sigma = (a * b < 0) ? -1: 1;
|
|
|
513 a = abs (a);
|
|
|
514 b = abs (b);
|
|
|
515 factor = sigma * (a / b);
|
|
|
516
|
|
|
517 lambda_matrix_row_add (S, n, i, i-1, -factor);
|
|
|
518 lambda_matrix_row_exchange (S, i, i-1);
|
|
|
519
|
|
|
520 lambda_matrix_row_add (U, m, i, i-1, -factor);
|
|
|
521 lambda_matrix_row_exchange (U, i, i-1);
|
|
|
522 }
|
|
|
523 }
|
|
|
524 }
|
|
|
525 }
|
|
|
526 }
|
|
|
527
|
|
|
528 /* Given an M x N integer matrix A, this function determines an M x M
|
|
|
529 unimodular matrix V, and an M x N echelon matrix S such that "A =
|
|
|
530 V.S". This decomposition is also known as "left Hermite".
|
|
|
531
|
|
|
532 Ref: Algorithm 2.2 page 36 in "Loop Transformations for
|
|
|
533 Restructuring Compilers" Utpal Banerjee. */
|
|
|
534
|
|
|
535 void
|
|
|
536 lambda_matrix_left_hermite (lambda_matrix A, int m, int n,
|
|
|
537 lambda_matrix S, lambda_matrix V)
|
|
|
538 {
|
|
|
539 int i, j, i0 = 0;
|
|
|
540
|
|
|
541 lambda_matrix_copy (A, S, m, n);
|
|
|
542 lambda_matrix_id (V, m);
|
|
|
543
|
|
|
544 for (j = 0; j < n; j++)
|
|
|
545 {
|
|
|
546 if (lambda_vector_first_nz (S[j], m, i0) < m)
|
|
|
547 {
|
|
|
548 ++i0;
|
|
|
549 for (i = m - 1; i >= i0; i--)
|
|
|
550 {
|
|
|
551 while (S[i][j] != 0)
|
|
|
552 {
|
|
|
553 int sigma, factor, a, b;
|
|
|
554
|
|
|
555 a = S[i-1][j];
|
|
|
556 b = S[i][j];
|
|
|
557 sigma = (a * b < 0) ? -1: 1;
|
|
|
558 a = abs (a);
|
|
|
559 b = abs (b);
|
|
|
560 factor = sigma * (a / b);
|
|
|
561
|
|
|
562 lambda_matrix_row_add (S, n, i, i-1, -factor);
|
|
|
563 lambda_matrix_row_exchange (S, i, i-1);
|
|
|
564
|
|
|
565 lambda_matrix_col_add (V, m, i-1, i, factor);
|
|
|
566 lambda_matrix_col_exchange (V, m, i, i-1);
|
|
|
567 }
|
|
|
568 }
|
|
|
569 }
|
|
|
570 }
|
|
|
571 }
|
|
|
572
|
|
|
573 /* When it exists, return the first nonzero row in MAT after row
|
|
|
574 STARTROW. Otherwise return rowsize. */
|
|
|
575
|
|
|
576 int
|
|
|
577 lambda_matrix_first_nz_vec (lambda_matrix mat, int rowsize, int colsize,
|
|
|
578 int startrow)
|
|
|
579 {
|
|
|
580 int j;
|
|
|
581 bool found = false;
|
|
|
582
|
|
|
583 for (j = startrow; (j < rowsize) && !found; j++)
|
|
|
584 {
|
|
|
585 if ((mat[j] != NULL)
|
|
|
586 && (lambda_vector_first_nz (mat[j], colsize, startrow) < colsize))
|
|
|
587 found = true;
|
|
|
588 }
|
|
|
589
|
|
|
590 if (found)
|
|
|
591 return j - 1;
|
|
|
592 return rowsize;
|
|
|
593 }
|
|
|
594
|
|
|
595 /* Calculate the projection of E sub k to the null space of B. */
|
|
|
596
|
|
|
597 void
|
|
|
598 lambda_matrix_project_to_null (lambda_matrix B, int rowsize,
|
|
|
599 int colsize, int k, lambda_vector x)
|
|
|
600 {
|
|
|
601 lambda_matrix M1, M2, M3, I;
|
|
|
602 int determinant;
|
|
|
603
|
|
|
604 /* Compute c(I-B^T inv(B B^T) B) e sub k. */
|
|
|
605
|
|
|
606 /* M1 is the transpose of B. */
|
|
|
607 M1 = lambda_matrix_new (colsize, colsize);
|
|
|
608 lambda_matrix_transpose (B, M1, rowsize, colsize);
|
|
|
609
|
|
|
610 /* M2 = B * B^T */
|
|
|
611 M2 = lambda_matrix_new (colsize, colsize);
|
|
|
612 lambda_matrix_mult (B, M1, M2, rowsize, colsize, rowsize);
|
|
|
613
|
|
|
614 /* M3 = inv(M2) */
|
|
|
615 M3 = lambda_matrix_new (colsize, colsize);
|
|
|
616 determinant = lambda_matrix_inverse (M2, M3, rowsize);
|
|
|
617
|
|
|
618 /* M2 = B^T (inv(B B^T)) */
|
|
|
619 lambda_matrix_mult (M1, M3, M2, colsize, rowsize, rowsize);
|
|
|
620
|
|
|
621 /* M1 = B^T (inv(B B^T)) B */
|
|
|
622 lambda_matrix_mult (M2, B, M1, colsize, rowsize, colsize);
|
|
|
623 lambda_matrix_negate (M1, M1, colsize, colsize);
|
|
|
624
|
|
|
625 I = lambda_matrix_new (colsize, colsize);
|
|
|
626 lambda_matrix_id (I, colsize);
|
|
|
627
|
|
|
628 lambda_matrix_add_mc (I, determinant, M1, 1, M2, colsize, colsize);
|
|
|
629
|
|
|
630 lambda_matrix_get_column (M2, colsize, k - 1, x);
|
|
|
631
|
|
|
632 }
|
|
|
633
|
|
|
634 /* Multiply a vector VEC by a matrix MAT.
|
|
|
635 MAT is an M*N matrix, and VEC is a vector with length N. The result
|
|
|
636 is stored in DEST which must be a vector of length M. */
|
|
|
637
|
|
|
638 void
|
|
|
639 lambda_matrix_vector_mult (lambda_matrix matrix, int m, int n,
|
|
|
640 lambda_vector vec, lambda_vector dest)
|
|
|
641 {
|
|
|
642 int i, j;
|
|
|
643
|
|
|
644 lambda_vector_clear (dest, m);
|
|
|
645 for (i = 0; i < m; i++)
|
|
|
646 for (j = 0; j < n; j++)
|
|
|
647 dest[i] += matrix[i][j] * vec[j];
|
|
|
648 }
|
|
|
649
|
|
|
650 /* Print out an M x N matrix MAT to OUTFILE. */
|
|
|
651
|
|
|
652 void
|
|
|
653 print_lambda_matrix (FILE * outfile, lambda_matrix matrix, int m, int n)
|
|
|
654 {
|
|
|
655 int i;
|
|
|
656
|
|
|
657 for (i = 0; i < m; i++)
|
|
|
658 print_lambda_vector (outfile, matrix[i], n);
|
|
|
659 fprintf (outfile, "\n");
|
|
|
660 }
|
|
|
661
|
