gcc/lambda-mat.c

   1 /* Integer matrix math routines
   2    Copyright (C) 2003, 2004, 2005 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 2, 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 COPYING.  If not, write to the Free
  19 Software Foundation, 51 Franklin Street, Fifth Floor, Boston, MA
  20 02110-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
  29 static void lambda_matrix_get_column (lambda_matrix, int, int,
  30                                       lambda_vector);
  31
  32 /* Allocate a matrix of M rows x  N cols.  */
  33
  34 lambda_matrix
  35 lambda_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
  50 void
  51 lambda_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
  62 void
  63 lambda_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
  74 bool
  75 lambda_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
  97 void
  98 lambda_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
 109 void
 110 lambda_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
 122 void
 123 lambda_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
 134 void
 135 lambda_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
 149 void
 150 lambda_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
 170 static void
 171 lambda_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
 182 void
 183 lambda_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
 198 void
 199 lambda_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
 211 void
 212 lambda_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
 225 void
 226 lambda_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
 233 void
 234 lambda_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
 244 void
 245 lambda_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
 260 void
 261 lambda_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
 274 void
 275 lambda_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
 285 void
 286 lambda_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
 309 static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int);
 310
 311 int
 312 lambda_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
 342 static int
 343 lambda_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
 441 void
 442 lambda_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
 489 void
 490 lambda_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
 534 void
 535 lambda_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
 575 int
 576 lambda_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
 596 void
 597 lambda_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
 637 void
 638 lambda_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
 651 void
 652 print_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