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