| /* Integer matrix math routines |
| Copyright (C) 2003, 2004 Free Software Foundation, Inc. |
| Contributed by Daniel Berlin <dberlin@dberlin.org>. |
| |
| This file is part of GCC. |
| |
| GCC is free software; you can redistribute it and/or modify it under |
| the terms of the GNU General Public License as published by the Free |
| Software Foundation; either version 2, or (at your option) any later |
| version. |
| |
| GCC is distributed in the hope that it will be useful, but WITHOUT ANY |
| WARRANTY; without even the implied warranty of MERCHANTABILITY or |
| FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License |
| for more details. |
| |
| You should have received a copy of the GNU General Public License |
| along with GCC; see the file COPYING. If not, write to the Free |
| Software Foundation, 59 Temple Place - Suite 330, Boston, MA |
| 02111-1307, USA. */ |
| #include "config.h" |
| #include "system.h" |
| #include "coretypes.h" |
| #include "tm.h" |
| #include "ggc.h" |
| #include "varray.h" |
| #include "tree.h" |
| #include "lambda.h" |
| |
| static void lambda_matrix_get_column (lambda_matrix, int, int, |
| lambda_vector); |
| |
| /* Allocate a matrix of M rows x N cols. */ |
| |
| lambda_matrix |
| lambda_matrix_new (int m, int n) |
| { |
| lambda_matrix mat; |
| int i; |
| |
| mat = ggc_alloc (m * sizeof (lambda_vector)); |
| |
| for (i = 0; i < m; i++) |
| mat[i] = lambda_vector_new (n); |
| |
| return mat; |
| } |
| |
| /* Copy the elements of M x N matrix MAT1 to MAT2. */ |
| |
| void |
| lambda_matrix_copy (lambda_matrix mat1, lambda_matrix mat2, |
| int m, int n) |
| { |
| int i; |
| |
| for (i = 0; i < m; i++) |
| lambda_vector_copy (mat1[i], mat2[i], n); |
| } |
| |
| /* Store the N x N identity matrix in MAT. */ |
| |
| void |
| lambda_matrix_id (lambda_matrix mat, int size) |
| { |
| int i, j; |
| |
| for (i = 0; i < size; i++) |
| for (j = 0; j < size; j++) |
| mat[i][j] = (i == j) ? 1 : 0; |
| } |
| |
| /* Return true if MAT is the identity matrix of SIZE */ |
| |
| bool |
| lambda_matrix_id_p (lambda_matrix mat, int size) |
| { |
| int i, j; |
| for (i = 0; i < size; i++) |
| for (j = 0; j < size; j++) |
| { |
| if (i == j) |
| { |
| if (mat[i][j] != 1) |
| return false; |
| } |
| else |
| { |
| if (mat[i][j] != 0) |
| return false; |
| } |
| } |
| return true; |
| } |
| |
| /* Negate the elements of the M x N matrix MAT1 and store it in MAT2. */ |
| |
| void |
| lambda_matrix_negate (lambda_matrix mat1, lambda_matrix mat2, int m, int n) |
| { |
| int i; |
| |
| for (i = 0; i < m; i++) |
| lambda_vector_negate (mat1[i], mat2[i], n); |
| } |
| |
| /* Take the transpose of matrix MAT1 and store it in MAT2. |
| MAT1 is an M x N matrix, so MAT2 must be N x M. */ |
| |
| void |
| lambda_matrix_transpose (lambda_matrix mat1, lambda_matrix mat2, int m, int n) |
| { |
| int i, j; |
| |
| for (i = 0; i < n; i++) |
| for (j = 0; j < m; j++) |
| mat2[i][j] = mat1[j][i]; |
| } |
| |
| |
| /* Add two M x N matrices together: MAT3 = MAT1+MAT2. */ |
| |
| void |
| lambda_matrix_add (lambda_matrix mat1, lambda_matrix mat2, |
| lambda_matrix mat3, int m, int n) |
| { |
| int i; |
| |
| for (i = 0; i < m; i++) |
| lambda_vector_add (mat1[i], mat2[i], mat3[i], n); |
| } |
| |
| /* MAT3 = CONST1 * MAT1 + CONST2 * MAT2. All matrices are M x N. */ |
| |
| void |
| lambda_matrix_add_mc (lambda_matrix mat1, int const1, |
| lambda_matrix mat2, int const2, |
| lambda_matrix mat3, int m, int n) |
| { |
| int i; |
| |
| for (i = 0; i < m; i++) |
| lambda_vector_add_mc (mat1[i], const1, mat2[i], const2, mat3[i], n); |
| } |
| |
| /* Multiply two matrices: MAT3 = MAT1 * MAT2. |
| MAT1 is an M x R matrix, and MAT2 is R x N. The resulting MAT2 |
| must therefore be M x N. */ |
| |
| void |
| lambda_matrix_mult (lambda_matrix mat1, lambda_matrix mat2, |
| lambda_matrix mat3, int m, int r, int n) |
| { |
| |
| int i, j, k; |
| |
| for (i = 0; i < m; i++) |
| { |
| for (j = 0; j < n; j++) |
| { |
| mat3[i][j] = 0; |
| for (k = 0; k < r; k++) |
| mat3[i][j] += mat1[i][k] * mat2[k][j]; |
| } |
| } |
| } |
| |
| /* Get column COL from the matrix MAT and store it in VEC. MAT has |
| N rows, so the length of VEC must be N. */ |
| |
| static void |
| lambda_matrix_get_column (lambda_matrix mat, int n, int col, |
| lambda_vector vec) |
| { |
| int i; |
| |
| for (i = 0; i < n; i++) |
| vec[i] = mat[i][col]; |
| } |
| |
| /* Delete rows r1 to r2 (not including r2). */ |
| |
| void |
| lambda_matrix_delete_rows (lambda_matrix mat, int rows, int from, int to) |
| { |
| int i; |
| int dist; |
| dist = to - from; |
| |
| for (i = to; i < rows; i++) |
| mat[i - dist] = mat[i]; |
| |
| for (i = rows - dist; i < rows; i++) |
| mat[i] = NULL; |
| } |
| |
| /* Swap rows R1 and R2 in matrix MAT. */ |
| |
| void |
| lambda_matrix_row_exchange (lambda_matrix mat, int r1, int r2) |
| { |
| lambda_vector row; |
| |
| row = mat[r1]; |
| mat[r1] = mat[r2]; |
| mat[r2] = row; |
| } |
| |
| /* Add a multiple of row R1 of matrix MAT with N columns to row R2: |
| R2 = R2 + CONST1 * R1. */ |
| |
| void |
| lambda_matrix_row_add (lambda_matrix mat, int n, int r1, int r2, int const1) |
| { |
| int i; |
| |
| if (const1 == 0) |
| return; |
| |
| for (i = 0; i < n; i++) |
| mat[r2][i] += const1 * mat[r1][i]; |
| } |
| |
| /* Negate row R1 of matrix MAT which has N columns. */ |
| |
| void |
| lambda_matrix_row_negate (lambda_matrix mat, int n, int r1) |
| { |
| lambda_vector_negate (mat[r1], mat[r1], n); |
| } |
| |
| /* Multiply row R1 of matrix MAT with N columns by CONST1. */ |
| |
| void |
| lambda_matrix_row_mc (lambda_matrix mat, int n, int r1, int const1) |
| { |
| int i; |
| |
| for (i = 0; i < n; i++) |
| mat[r1][i] *= const1; |
| } |
| |
| /* Exchange COL1 and COL2 in matrix MAT. M is the number of rows. */ |
| |
| void |
| lambda_matrix_col_exchange (lambda_matrix mat, int m, int col1, int col2) |
| { |
| int i; |
| int tmp; |
| for (i = 0; i < m; i++) |
| { |
| tmp = mat[i][col1]; |
| mat[i][col1] = mat[i][col2]; |
| mat[i][col2] = tmp; |
| } |
| } |
| |
| /* Add a multiple of column C1 of matrix MAT with M rows to column C2: |
| C2 = C2 + CONST1 * C1. */ |
| |
| void |
| lambda_matrix_col_add (lambda_matrix mat, int m, int c1, int c2, int const1) |
| { |
| int i; |
| |
| if (const1 == 0) |
| return; |
| |
| for (i = 0; i < m; i++) |
| mat[i][c2] += const1 * mat[i][c1]; |
| } |
| |
| /* Negate column C1 of matrix MAT which has M rows. */ |
| |
| void |
| lambda_matrix_col_negate (lambda_matrix mat, int m, int c1) |
| { |
| int i; |
| |
| for (i = 0; i < m; i++) |
| mat[i][c1] *= -1; |
| } |
| |
| /* Multiply column C1 of matrix MAT with M rows by CONST1. */ |
| |
| void |
| lambda_matrix_col_mc (lambda_matrix mat, int m, int c1, int const1) |
| { |
| int i; |
| |
| for (i = 0; i < m; i++) |
| mat[i][c1] *= const1; |
| } |
| |
| /* Compute the inverse of the N x N matrix MAT and store it in INV. |
| |
| We don't _really_ compute the inverse of MAT. Instead we compute |
| det(MAT)*inv(MAT), and we return det(MAT) to the caller as the function |
| result. This is necessary to preserve accuracy, because we are dealing |
| with integer matrices here. |
| |
| The algorithm used here is a column based Gauss-Jordan elimination on MAT |
| and the identity matrix in parallel. The inverse is the result of applying |
| the same operations on the identity matrix that reduce MAT to the identity |
| matrix. |
| |
| When MAT is a 2 x 2 matrix, we don't go through the whole process, because |
| it is easily inverted by inspection and it is a very common case. */ |
| |
| static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int); |
| |
| int |
| lambda_matrix_inverse (lambda_matrix mat, lambda_matrix inv, int n) |
| { |
| if (n == 2) |
| { |
| int a, b, c, d, det; |
| a = mat[0][0]; |
| b = mat[1][0]; |
| c = mat[0][1]; |
| d = mat[1][1]; |
| inv[0][0] = d; |
| inv[0][1] = -c; |
| inv[1][0] = -b; |
| inv[1][1] = a; |
| det = (a * d - b * c); |
| if (det < 0) |
| { |
| det *= -1; |
| inv[0][0] *= -1; |
| inv[1][0] *= -1; |
| inv[0][1] *= -1; |
| inv[1][1] *= -1; |
| } |
| return det; |
| } |
| else |
| return lambda_matrix_inverse_hard (mat, inv, n); |
| } |
| |
| /* If MAT is not a special case, invert it the hard way. */ |
| |
| static int |
| lambda_matrix_inverse_hard (lambda_matrix mat, lambda_matrix inv, int n) |
| { |
| lambda_vector row; |
| lambda_matrix temp; |
| int i, j; |
| int determinant; |
| |
| temp = lambda_matrix_new (n, n); |
| lambda_matrix_copy (mat, temp, n, n); |
| lambda_matrix_id (inv, n); |
| |
| /* Reduce TEMP to a lower triangular form, applying the same operations on |
| INV which starts as the identity matrix. N is the number of rows and |
| columns. */ |
| for (j = 0; j < n; j++) |
| { |
| row = temp[j]; |
| |
| /* Make every element in the current row positive. */ |
| for (i = j; i < n; i++) |
| if (row[i] < 0) |
| { |
| lambda_matrix_col_negate (temp, n, i); |
| lambda_matrix_col_negate (inv, n, i); |
| } |
| |
| /* Sweep the upper triangle. Stop when only the diagonal element in the |
| current row is nonzero. */ |
| while (lambda_vector_first_nz (row, n, j + 1) < n) |
| { |
| int min_col = lambda_vector_min_nz (row, n, j); |
| lambda_matrix_col_exchange (temp, n, j, min_col); |
| lambda_matrix_col_exchange (inv, n, j, min_col); |
| |
| for (i = j + 1; i < n; i++) |
| { |
| int factor; |
| |
| factor = -1 * row[i]; |
| if (row[j] != 1) |
| factor /= row[j]; |
| |
| lambda_matrix_col_add (temp, n, j, i, factor); |
| lambda_matrix_col_add (inv, n, j, i, factor); |
| } |
| } |
| } |
| |
| /* Reduce TEMP from a lower triangular to the identity matrix. Also compute |
| the determinant, which now is simply the product of the elements on the |
| diagonal of TEMP. If one of these elements is 0, the matrix has 0 as an |
| eigenvalue so it is singular and hence not invertible. */ |
| determinant = 1; |
| for (j = n - 1; j >= 0; j--) |
| { |
| int diagonal; |
| |
| row = temp[j]; |
| diagonal = row[j]; |
| |
| /* If the matrix is singular, abort. */ |
| if (diagonal == 0) |
| abort (); |
| |
| determinant = determinant * diagonal; |
| |
| /* If the diagonal is not 1, then multiply the each row by the |
| diagonal so that the middle number is now 1, rather than a |
| rational. */ |
| if (diagonal != 1) |
| { |
| for (i = 0; i < j; i++) |
| lambda_matrix_col_mc (inv, n, i, diagonal); |
| for (i = j + 1; i < n; i++) |
| lambda_matrix_col_mc (inv, n, i, diagonal); |
| |
| row[j] = diagonal = 1; |
| } |
| |
| /* Sweep the lower triangle column wise. */ |
| for (i = j - 1; i >= 0; i--) |
| { |
| if (row[i]) |
| { |
| int factor = -row[i]; |
| lambda_matrix_col_add (temp, n, j, i, factor); |
| lambda_matrix_col_add (inv, n, j, i, factor); |
| } |
| |
| } |
| } |
| |
| return determinant; |
| } |
| |
| /* Decompose a N x N matrix MAT to a product of a lower triangular H |
| and a unimodular U matrix such that MAT = H.U. N is the size of |
| the rows of MAT. */ |
| |
| void |
| lambda_matrix_hermite (lambda_matrix mat, int n, |
| lambda_matrix H, lambda_matrix U) |
| { |
| lambda_vector row; |
| int i, j, factor, minimum_col; |
| |
| lambda_matrix_copy (mat, H, n, n); |
| lambda_matrix_id (U, n); |
| |
| for (j = 0; j < n; j++) |
| { |
| row = H[j]; |
| |
| /* Make every element of H[j][j..n] positive. */ |
| for (i = j; i < n; i++) |
| { |
| if (row[i] < 0) |
| { |
| lambda_matrix_col_negate (H, n, i); |
| lambda_vector_negate (U[i], U[i], n); |
| } |
| } |
| |
| /* Stop when only the diagonal element is nonzero. */ |
| while (lambda_vector_first_nz (row, n, j + 1) < n) |
| { |
| minimum_col = lambda_vector_min_nz (row, n, j); |
| lambda_matrix_col_exchange (H, n, j, minimum_col); |
| lambda_matrix_row_exchange (U, j, minimum_col); |
| |
| for (i = j + 1; i < n; i++) |
| { |
| factor = row[i] / row[j]; |
| lambda_matrix_col_add (H, n, j, i, -1 * factor); |
| lambda_matrix_row_add (U, n, i, j, factor); |
| } |
| } |
| } |
| } |
| |
| /* Given an M x N integer matrix A, this function determines an M x |
| M unimodular matrix U, and an M x N echelon matrix S such that |
| "U.A = S". This decomposition is also known as "right Hermite". |
| |
| Ref: Algorithm 2.1 page 33 in "Loop Transformations for |
| Restructuring Compilers" Utpal Banerjee. */ |
| |
| void |
| lambda_matrix_right_hermite (lambda_matrix A, int m, int n, |
| lambda_matrix S, lambda_matrix U) |
| { |
| int i, j, i0 = 0; |
| |
| lambda_matrix_copy (A, S, m, n); |
| lambda_matrix_id (U, m); |
| |
| for (j = 0; j < n; j++) |
| { |
| if (lambda_vector_first_nz (S[j], m, i0) < m) |
| { |
| ++i0; |
| for (i = m - 1; i >= i0; i--) |
| { |
| while (S[i][j] != 0) |
| { |
| int sigma, factor, a, b; |
| |
| a = S[i-1][j]; |
| b = S[i][j]; |
| sigma = (a * b < 0) ? -1: 1; |
| a = abs (a); |
| b = abs (b); |
| factor = sigma * (a / b); |
| |
| lambda_matrix_row_add (S, n, i, i-1, -factor); |
| lambda_matrix_row_exchange (S, i, i-1); |
| |
| lambda_matrix_row_add (U, m, i, i-1, -factor); |
| lambda_matrix_row_exchange (U, i, i-1); |
| } |
| } |
| } |
| } |
| } |
| |
| /* Given an M x N integer matrix A, this function determines an M x M |
| unimodular matrix V, and an M x N echelon matrix S such that "A = |
| V.S". This decomposition is also known as "left Hermite". |
| |
| Ref: Algorithm 2.2 page 36 in "Loop Transformations for |
| Restructuring Compilers" Utpal Banerjee. */ |
| |
| void |
| lambda_matrix_left_hermite (lambda_matrix A, int m, int n, |
| lambda_matrix S, lambda_matrix V) |
| { |
| int i, j, i0 = 0; |
| |
| lambda_matrix_copy (A, S, m, n); |
| lambda_matrix_id (V, m); |
| |
| for (j = 0; j < n; j++) |
| { |
| if (lambda_vector_first_nz (S[j], m, i0) < m) |
| { |
| ++i0; |
| for (i = m - 1; i >= i0; i--) |
| { |
| while (S[i][j] != 0) |
| { |
| int sigma, factor, a, b; |
| |
| a = S[i-1][j]; |
| b = S[i][j]; |
| sigma = (a * b < 0) ? -1: 1; |
| a = abs (a); |
| b = abs (b); |
| factor = sigma * (a / b); |
| |
| lambda_matrix_row_add (S, n, i, i-1, -factor); |
| lambda_matrix_row_exchange (S, i, i-1); |
| |
| lambda_matrix_col_add (V, m, i-1, i, factor); |
| lambda_matrix_col_exchange (V, m, i, i-1); |
| } |
| } |
| } |
| } |
| } |
| |
| /* When it exists, return the first nonzero row in MAT after row |
| STARTROW. Otherwise return rowsize. */ |
| |
| int |
| lambda_matrix_first_nz_vec (lambda_matrix mat, int rowsize, int colsize, |
| int startrow) |
| { |
| int j; |
| bool found = false; |
| |
| for (j = startrow; (j < rowsize) && !found; j++) |
| { |
| if ((mat[j] != NULL) |
| && (lambda_vector_first_nz (mat[j], colsize, startrow) < colsize)) |
| found = true; |
| } |
| |
| if (found) |
| return j - 1; |
| return rowsize; |
| } |
| |
| /* Calculate the projection of E sub k to the null space of B. */ |
| |
| void |
| lambda_matrix_project_to_null (lambda_matrix B, int rowsize, |
| int colsize, int k, lambda_vector x) |
| { |
| lambda_matrix M1, M2, M3, I; |
| int determinant; |
| |
| /* Compute c(I-B^T inv(B B^T) B) e sub k. */ |
| |
| /* M1 is the transpose of B. */ |
| M1 = lambda_matrix_new (colsize, colsize); |
| lambda_matrix_transpose (B, M1, rowsize, colsize); |
| |
| /* M2 = B * B^T */ |
| M2 = lambda_matrix_new (colsize, colsize); |
| lambda_matrix_mult (B, M1, M2, rowsize, colsize, rowsize); |
| |
| /* M3 = inv(M2) */ |
| M3 = lambda_matrix_new (colsize, colsize); |
| determinant = lambda_matrix_inverse (M2, M3, rowsize); |
| |
| /* M2 = B^T (inv(B B^T)) */ |
| lambda_matrix_mult (M1, M3, M2, colsize, rowsize, rowsize); |
| |
| /* M1 = B^T (inv(B B^T)) B */ |
| lambda_matrix_mult (M2, B, M1, colsize, rowsize, colsize); |
| lambda_matrix_negate (M1, M1, colsize, colsize); |
| |
| I = lambda_matrix_new (colsize, colsize); |
| lambda_matrix_id (I, colsize); |
| |
| lambda_matrix_add_mc (I, determinant, M1, 1, M2, colsize, colsize); |
| |
| lambda_matrix_get_column (M2, colsize, k - 1, x); |
| |
| } |
| |
| /* Multiply a vector VEC by a matrix MAT. |
| MAT is an M*N matrix, and VEC is a vector with length N. The result |
| is stored in DEST which must be a vector of length M. */ |
| |
| void |
| lambda_matrix_vector_mult (lambda_matrix matrix, int m, int n, |
| lambda_vector vec, lambda_vector dest) |
| { |
| int i, j; |
| |
| lambda_vector_clear (dest, m); |
| for (i = 0; i < m; i++) |
| for (j = 0; j < n; j++) |
| dest[i] += matrix[i][j] * vec[j]; |
| } |
| |
| /* Print out an M x N matrix MAT to OUTFILE. */ |
| |
| void |
| print_lambda_matrix (FILE * outfile, lambda_matrix matrix, int m, int n) |
| { |
| int i; |
| |
| for (i = 0; i < m; i++) |
| print_lambda_vector (outfile, matrix[i], n); |
| fprintf (outfile, "\n"); |
| } |
| |