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