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 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 3, 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 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 30static void lambda_matrix_get_column (lambda_matrix, int, int, 31 lambda_vector); 32 33/* Allocate a matrix of M rows x N cols. */ 34 35lambda_matrix 36lambda_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 51void 52lambda_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 63void 64lambda_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 75bool 76lambda_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 98void 99lambda_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 110void 111lambda_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 123void 124lambda_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 135void 136lambda_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 150void 151lambda_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 171static void 172lambda_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 183void 184lambda_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 199void 200lambda_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 212void 213lambda_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 226void 227lambda_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 234void 235lambda_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 245void 246lambda_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 261void 262lambda_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 275void 276lambda_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 286void 287lambda_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 310static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int); 311 312int 313lambda_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 343static int 344lambda_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 442void 443lambda_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 490void 491lambda_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 535void 536lambda_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 576int 577lambda_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 597void 598lambda_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 638void 639lambda_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 652void 653print_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 662