gcc/lambda-mat.c

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