gcc/lambda-mat.c

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