1 /* Integer matrix math routines
2 Copyright (C) 2003, 2004, 2005, 2007, 2008, 2010
3 Free Software Foundation, Inc.
4 Contributed by Daniel Berlin <dberlin@dberlin.org>.
6 This file is part of GCC.
8 GCC is free software; you can redistribute it and/or modify it under
9 the terms of the GNU General Public License as published by the Free
10 Software Foundation; either version 3, or (at your option) any later
13 GCC is distributed in the hope that it will be useful, but WITHOUT ANY
14 WARRANTY; without even the implied warranty of MERCHANTABILITY or
15 FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License
18 You should have received a copy of the GNU General Public License
19 along with GCC; see the file COPYING3. If not see
20 <http://www.gnu.org/licenses/>. */
24 #include "coretypes.h"
25 #include "tree-flow.h"
28 /* Allocate a matrix of M rows x N cols. */
31 lambda_matrix_new (int m, int n, struct obstack * lambda_obstack)
36 mat = (lambda_matrix) obstack_alloc (lambda_obstack,
37 sizeof (lambda_vector *) * m);
39 for (i = 0; i < m; i++)
40 mat[i] = lambda_vector_new (n);
45 /* Copy the elements of M x N matrix MAT1 to MAT2. */
48 lambda_matrix_copy (lambda_matrix mat1, lambda_matrix mat2,
53 for (i = 0; i < m; i++)
54 lambda_vector_copy (mat1[i], mat2[i], n);
57 /* Store the N x N identity matrix in MAT. */
60 lambda_matrix_id (lambda_matrix mat, int size)
64 for (i = 0; i < size; i++)
65 for (j = 0; j < size; j++)
66 mat[i][j] = (i == j) ? 1 : 0;
69 /* Return true if MAT is the identity matrix of SIZE */
72 lambda_matrix_id_p (lambda_matrix mat, int size)
75 for (i = 0; i < size; i++)
76 for (j = 0; j < size; j++)
92 /* Negate the elements of the M x N matrix MAT1 and store it in MAT2. */
95 lambda_matrix_negate (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
99 for (i = 0; i < m; i++)
100 lambda_vector_negate (mat1[i], mat2[i], n);
103 /* Take the transpose of matrix MAT1 and store it in MAT2.
104 MAT1 is an M x N matrix, so MAT2 must be N x M. */
107 lambda_matrix_transpose (lambda_matrix mat1, lambda_matrix mat2, int m, int n)
111 for (i = 0; i < n; i++)
112 for (j = 0; j < m; j++)
113 mat2[i][j] = mat1[j][i];
117 /* Add two M x N matrices together: MAT3 = MAT1+MAT2. */
120 lambda_matrix_add (lambda_matrix mat1, lambda_matrix mat2,
121 lambda_matrix mat3, int m, int n)
125 for (i = 0; i < m; i++)
126 lambda_vector_add (mat1[i], mat2[i], mat3[i], n);
129 /* MAT3 = CONST1 * MAT1 + CONST2 * MAT2. All matrices are M x N. */
132 lambda_matrix_add_mc (lambda_matrix mat1, int const1,
133 lambda_matrix mat2, int const2,
134 lambda_matrix mat3, int m, int n)
138 for (i = 0; i < m; i++)
139 lambda_vector_add_mc (mat1[i], const1, mat2[i], const2, mat3[i], n);
142 /* Multiply two matrices: MAT3 = MAT1 * MAT2.
143 MAT1 is an M x R matrix, and MAT2 is R x N. The resulting MAT2
144 must therefore be M x N. */
147 lambda_matrix_mult (lambda_matrix mat1, lambda_matrix mat2,
148 lambda_matrix mat3, int m, int r, int n)
153 for (i = 0; i < m; i++)
155 for (j = 0; j < n; j++)
158 for (k = 0; k < r; k++)
159 mat3[i][j] += mat1[i][k] * mat2[k][j];
164 /* Delete rows r1 to r2 (not including r2). */
167 lambda_matrix_delete_rows (lambda_matrix mat, int rows, int from, int to)
173 for (i = to; i < rows; i++)
174 mat[i - dist] = mat[i];
176 for (i = rows - dist; i < rows; i++)
180 /* Swap rows R1 and R2 in matrix MAT. */
183 lambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2)
192 /* Add a multiple of row R1 of matrix MAT with N columns to row R2:
193 R2 = R2 + CONST1 * R1. */
196 lambda_matrix_row_add (lambda_matrix mat, int n, int r1, int r2, int const1)
203 for (i = 0; i < n; i++)
204 mat[r2][i] += const1 * mat[r1][i];
207 /* Negate row R1 of matrix MAT which has N columns. */
210 lambda_matrix_row_negate (lambda_matrix mat, int n, int r1)
212 lambda_vector_negate (mat[r1], mat[r1], n);
215 /* Multiply row R1 of matrix MAT with N columns by CONST1. */
218 lambda_matrix_row_mc (lambda_matrix mat, int n, int r1, int const1)
222 for (i = 0; i < n; i++)
223 mat[r1][i] *= const1;
226 /* Exchange COL1 and COL2 in matrix MAT. M is the number of rows. */
229 lambda_matrix_col_exchange (lambda_matrix mat, int m, int col1, int col2)
233 for (i = 0; i < m; i++)
236 mat[i][col1] = mat[i][col2];
241 /* Add a multiple of column C1 of matrix MAT with M rows to column C2:
242 C2 = C2 + CONST1 * C1. */
245 lambda_matrix_col_add (lambda_matrix mat, int m, int c1, int c2, int const1)
252 for (i = 0; i < m; i++)
253 mat[i][c2] += const1 * mat[i][c1];
256 /* Negate column C1 of matrix MAT which has M rows. */
259 lambda_matrix_col_negate (lambda_matrix mat, int m, int c1)
263 for (i = 0; i < m; i++)
267 /* Multiply column C1 of matrix MAT with M rows by CONST1. */
270 lambda_matrix_col_mc (lambda_matrix mat, int m, int c1, int const1)
274 for (i = 0; i < m; i++)
275 mat[i][c1] *= const1;
278 /* Compute the inverse of the N x N matrix MAT and store it in INV.
280 We don't _really_ compute the inverse of MAT. Instead we compute
281 det(MAT)*inv(MAT), and we return det(MAT) to the caller as the function
282 result. This is necessary to preserve accuracy, because we are dealing
283 with integer matrices here.
285 The algorithm used here is a column based Gauss-Jordan elimination on MAT
286 and the identity matrix in parallel. The inverse is the result of applying
287 the same operations on the identity matrix that reduce MAT to the identity
290 When MAT is a 2 x 2 matrix, we don't go through the whole process, because
291 it is easily inverted by inspection and it is a very common case. */
293 static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int,
297 lambda_matrix_inverse (lambda_matrix mat, lambda_matrix inv, int n,
298 struct obstack * lambda_obstack)
311 det = (a * d - b * c);
323 return lambda_matrix_inverse_hard (mat, inv, n, lambda_obstack);
326 /* If MAT is not a special case, invert it the hard way. */
329 lambda_matrix_inverse_hard (lambda_matrix mat, lambda_matrix inv, int n,
330 struct obstack * lambda_obstack)
337 temp = lambda_matrix_new (n, n, lambda_obstack);
338 lambda_matrix_copy (mat, temp, n, n);
339 lambda_matrix_id (inv, n);
341 /* Reduce TEMP to a lower triangular form, applying the same operations on
342 INV which starts as the identity matrix. N is the number of rows and
344 for (j = 0; j < n; j++)
348 /* Make every element in the current row positive. */
349 for (i = j; i < n; i++)
352 lambda_matrix_col_negate (temp, n, i);
353 lambda_matrix_col_negate (inv, n, i);
356 /* Sweep the upper triangle. Stop when only the diagonal element in the
357 current row is nonzero. */
358 while (lambda_vector_first_nz (row, n, j + 1) < n)
360 int min_col = lambda_vector_min_nz (row, n, j);
361 lambda_matrix_col_exchange (temp, n, j, min_col);
362 lambda_matrix_col_exchange (inv, n, j, min_col);
364 for (i = j + 1; i < n; i++)
368 factor = -1 * row[i];
372 lambda_matrix_col_add (temp, n, j, i, factor);
373 lambda_matrix_col_add (inv, n, j, i, factor);
378 /* Reduce TEMP from a lower triangular to the identity matrix. Also compute
379 the determinant, which now is simply the product of the elements on the
380 diagonal of TEMP. If one of these elements is 0, the matrix has 0 as an
381 eigenvalue so it is singular and hence not invertible. */
383 for (j = n - 1; j >= 0; j--)
390 /* The matrix must not be singular. */
391 gcc_assert (diagonal);
393 determinant = determinant * diagonal;
395 /* If the diagonal is not 1, then multiply the each row by the
396 diagonal so that the middle number is now 1, rather than a
400 for (i = 0; i < j; i++)
401 lambda_matrix_col_mc (inv, n, i, diagonal);
402 for (i = j + 1; i < n; i++)
403 lambda_matrix_col_mc (inv, n, i, diagonal);
405 row[j] = diagonal = 1;
408 /* Sweep the lower triangle column wise. */
409 for (i = j - 1; i >= 0; i--)
413 int factor = -row[i];
414 lambda_matrix_col_add (temp, n, j, i, factor);
415 lambda_matrix_col_add (inv, n, j, i, factor);
424 /* Decompose a N x N matrix MAT to a product of a lower triangular H
425 and a unimodular U matrix such that MAT = H.U. N is the size of
429 lambda_matrix_hermite (lambda_matrix mat, int n,
430 lambda_matrix H, lambda_matrix U)
433 int i, j, factor, minimum_col;
435 lambda_matrix_copy (mat, H, n, n);
436 lambda_matrix_id (U, n);
438 for (j = 0; j < n; j++)
442 /* Make every element of H[j][j..n] positive. */
443 for (i = j; i < n; i++)
447 lambda_matrix_col_negate (H, n, i);
448 lambda_vector_negate (U[i], U[i], n);
452 /* Stop when only the diagonal element is nonzero. */
453 while (lambda_vector_first_nz (row, n, j + 1) < n)
455 minimum_col = lambda_vector_min_nz (row, n, j);
456 lambda_matrix_col_exchange (H, n, j, minimum_col);
457 lambda_matrix_row_exchange (U, j, minimum_col);
459 for (i = j + 1; i < n; i++)
461 factor = row[i] / row[j];
462 lambda_matrix_col_add (H, n, j, i, -1 * factor);
463 lambda_matrix_row_add (U, n, i, j, factor);
469 /* Given an M x N integer matrix A, this function determines an M x
470 M unimodular matrix U, and an M x N echelon matrix S such that
471 "U.A = S". This decomposition is also known as "right Hermite".
473 Ref: Algorithm 2.1 page 33 in "Loop Transformations for
474 Restructuring Compilers" Utpal Banerjee. */
477 lambda_matrix_right_hermite (lambda_matrix A, int m, int n,
478 lambda_matrix S, lambda_matrix U)
482 lambda_matrix_copy (A, S, m, n);
483 lambda_matrix_id (U, m);
485 for (j = 0; j < n; j++)
487 if (lambda_vector_first_nz (S[j], m, i0) < m)
490 for (i = m - 1; i >= i0; i--)
494 int sigma, factor, a, b;
498 sigma = (a * b < 0) ? -1: 1;
501 factor = sigma * (a / b);
503 lambda_matrix_row_add (S, n, i, i-1, -factor);
504 lambda_matrix_row_exchange (S, i, i-1);
506 lambda_matrix_row_add (U, m, i, i-1, -factor);
507 lambda_matrix_row_exchange (U, i, i-1);
514 /* Given an M x N integer matrix A, this function determines an M x M
515 unimodular matrix V, and an M x N echelon matrix S such that "A =
516 V.S". This decomposition is also known as "left Hermite".
518 Ref: Algorithm 2.2 page 36 in "Loop Transformations for
519 Restructuring Compilers" Utpal Banerjee. */
522 lambda_matrix_left_hermite (lambda_matrix A, int m, int n,
523 lambda_matrix S, lambda_matrix V)
527 lambda_matrix_copy (A, S, m, n);
528 lambda_matrix_id (V, m);
530 for (j = 0; j < n; j++)
532 if (lambda_vector_first_nz (S[j], m, i0) < m)
535 for (i = m - 1; i >= i0; i--)
539 int sigma, factor, a, b;
543 sigma = (a * b < 0) ? -1: 1;
546 factor = sigma * (a / b);
548 lambda_matrix_row_add (S, n, i, i-1, -factor);
549 lambda_matrix_row_exchange (S, i, i-1);
551 lambda_matrix_col_add (V, m, i-1, i, factor);
552 lambda_matrix_col_exchange (V, m, i, i-1);
559 /* When it exists, return the first nonzero row in MAT after row
560 STARTROW. Otherwise return rowsize. */
563 lambda_matrix_first_nz_vec (lambda_matrix mat, int rowsize, int colsize,
569 for (j = startrow; (j < rowsize) && !found; j++)
572 && (lambda_vector_first_nz (mat[j], colsize, startrow) < colsize))
581 /* Multiply a vector VEC by a matrix MAT.
582 MAT is an M*N matrix, and VEC is a vector with length N. The result
583 is stored in DEST which must be a vector of length M. */
586 lambda_matrix_vector_mult (lambda_matrix matrix, int m, int n,
587 lambda_vector vec, lambda_vector dest)
591 lambda_vector_clear (dest, m);
592 for (i = 0; i < m; i++)
593 for (j = 0; j < n; j++)
594 dest[i] += matrix[i][j] * vec[j];
597 /* Print out an M x N matrix MAT to OUTFILE. */
600 print_lambda_matrix (FILE * outfile, lambda_matrix matrix, int m, int n)
604 for (i = 0; i < m; i++)
605 print_lambda_vector (outfile, matrix[i], n);
606 fprintf (outfile, "\n");