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 /* Return true if MAT is the identity matrix of SIZE */
  74
  75 bool
  76 lambda_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
  98 void
  99 lambda_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
 110 void
 111 lambda_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
 123 void
 124 lambda_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
 135 void
 136 lambda_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
 150 void
 151 lambda_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
 171 static void
 172 lambda_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
 183 void
 184 lambda_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
 199 void
 200 lambda_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
 212 void
 213 lambda_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
 226 void
 227 lambda_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
 234 void
 235 lambda_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
 245 void
 246 lambda_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
 261 void
 262 lambda_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
 275 void
 276 lambda_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
 286 void
 287 lambda_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
 310 static int lambda_matrix_inverse_hard (lambda_matrix, lambda_matrix, int);
 311
 312 int
 313 lambda_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
 343 static int
 344 lambda_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       /* If the matrix is singular, abort.  */
 405       if (diagonal == 0)
 406         abort ();
 407
 408       determinant = determinant * diagonal;
 409
 410       /* If the diagonal is not 1, then multiply the each row by the
 411          diagonal so that the middle number is now 1, rather than a
 412          rational.  */
 413       if (diagonal != 1)
 414         {
 415           for (i = 0; i < j; i++)
 416             lambda_matrix_col_mc (inv, n, i, diagonal);
 417           for (i = j + 1; i < n; i++)
 418             lambda_matrix_col_mc (inv, n, i, diagonal);
 419
 420           row[j] = diagonal = 1;
 421         }
 422
 423       /* Sweep the lower triangle column wise.  */
 424       for (i = j - 1; i >= 0; i--)
 425         {
 426           if (row[i])
 427             {
 428               int factor = -row[i];
 429               lambda_matrix_col_add (temp, n, j, i, factor);
 430               lambda_matrix_col_add (inv, n, j, i, factor);
 431             }
 432
 433         }
 434     }
 435
 436   return determinant;
 437 }
 438
 439 /* Decompose a N x N matrix MAT to a product of a lower triangular H
 440    and a unimodular U matrix such that MAT = H.U.  N is the size of
 441    the rows of MAT.  */
 442
 443 void
 444 lambda_matrix_hermite (lambda_matrix mat, int n,
 445                        lambda_matrix H, lambda_matrix U)
 446 {
 447   lambda_vector row;
 448   int i, j, factor, minimum_col;
 449
 450   lambda_matrix_copy (mat, H, n, n);
 451   lambda_matrix_id (U, n);
 452
 453   for (j = 0; j < n; j++)
 454     {
 455       row = H[j];
 456
 457       /* Make every element of H[j][j..n] positive.  */
 458       for (i = j; i < n; i++)
 459         {
 460           if (row[i] < 0)
 461             {
 462               lambda_matrix_col_negate (H, n, i);
 463               lambda_vector_negate (U[i], U[i], n);
 464             }
 465         }
 466
 467       /* Stop when only the diagonal element is nonzero.  */
 468       while (lambda_vector_first_nz (row, n, j + 1) < n)
 469         {
 470           minimum_col = lambda_vector_min_nz (row, n, j);
 471           lambda_matrix_col_exchange (H, n, j, minimum_col);
 472           lambda_matrix_row_exchange (U, j, minimum_col);
 473
 474           for (i = j + 1; i < n; i++)
 475             {
 476               factor = row[i] / row[j];
 477               lambda_matrix_col_add (H, n, j, i, -1 * factor);
 478               lambda_matrix_row_add (U, n, i, j, factor);
 479             }
 480         }
 481     }
 482 }
 483
 484 /* Given an M x N integer matrix A, this function determines an M x
 485    M unimodular matrix U, and an M x N echelon matrix S such that
 486    "U.A = S".  This decomposition is also known as "right Hermite".
 487
 488    Ref: Algorithm 2.1 page 33 in "Loop Transformations for
 489    Restructuring Compilers" Utpal Banerjee.  */
 490
 491 void
 492 lambda_matrix_right_hermite (lambda_matrix A, int m, int n,
 493                              lambda_matrix S, lambda_matrix U)
 494 {
 495   int i, j, i0 = 0;
 496
 497   lambda_matrix_copy (A, S, m, n);
 498   lambda_matrix_id (U, m);
 499
 500   for (j = 0; j < n; j++)
 501     {
 502       if (lambda_vector_first_nz (S[j], m, i0) < m)
 503         {
 504           ++i0;
 505           for (i = m - 1; i >= i0; i--)
 506             {
 507               while (S[i][j] != 0)
 508                 {
 509                   int sigma, factor, a, b;
 510
 511                   a = S[i-1][j];
 512                   b = S[i][j];
 513                   sigma = (a * b < 0) ? -1: 1;
 514                   a = abs (a);
 515                   b = abs (b);
 516                   factor = sigma * (a / b);
 517
 518                   lambda_matrix_row_add (S, n, i, i-1, -factor);
 519                   lambda_matrix_row_exchange (S, i, i-1);
 520
 521                   lambda_matrix_row_add (U, m, i, i-1, -factor);
 522                   lambda_matrix_row_exchange (U, i, i-1);
 523                 }
 524             }
 525         }
 526     }
 527 }
 528
 529 /* Given an M x N integer matrix A, this function determines an M x M
 530    unimodular matrix V, and an M x N echelon matrix S such that "A =
 531    V.S".  This decomposition is also known as "left Hermite".
 532
 533    Ref: Algorithm 2.2 page 36 in "Loop Transformations for
 534    Restructuring Compilers" Utpal Banerjee.  */
 535
 536 void
 537 lambda_matrix_left_hermite (lambda_matrix A, int m, int n,
 538                              lambda_matrix S, lambda_matrix V)
 539 {
 540   int i, j, i0 = 0;
 541
 542   lambda_matrix_copy (A, S, m, n);
 543   lambda_matrix_id (V, m);
 544
 545   for (j = 0; j < n; j++)
 546     {
 547       if (lambda_vector_first_nz (S[j], m, i0) < m)
 548         {
 549           ++i0;
 550           for (i = m - 1; i >= i0; i--)
 551             {
 552               while (S[i][j] != 0)
 553                 {
 554                   int sigma, factor, a, b;
 555
 556                   a = S[i-1][j];
 557                   b = S[i][j];
 558                   sigma = (a * b < 0) ? -1: 1;
 559                   a = abs (a);
 560       b = abs (b);
 561                   factor = sigma * (a / b);
 562
 563                   lambda_matrix_row_add (S, n, i, i-1, -factor);
 564                   lambda_matrix_row_exchange (S, i, i-1);
 565
 566                   lambda_matrix_col_add (V, m, i-1, i, factor);
 567                   lambda_matrix_col_exchange (V, m, i, i-1);
 568                 }
 569             }
 570         }
 571     }
 572 }
 573
 574 /* When it exists, return the first nonzero row in MAT after row
 575    STARTROW.  Otherwise return rowsize.  */
 576
 577 int
 578 lambda_matrix_first_nz_vec (lambda_matrix mat, int rowsize, int colsize,
 579                             int startrow)
 580 {
 581   int j;
 582   bool found = false;
 583
 584   for (j = startrow; (j < rowsize) && !found; j++)
 585     {
 586       if ((mat[j] != NULL)
 587           && (lambda_vector_first_nz (mat[j], colsize, startrow) < colsize))
 588         found = true;
 589     }
 590
 591   if (found)
 592     return j - 1;
 593   return rowsize;
 594 }
 595
 596 /* Calculate the projection of E sub k to the null space of B.  */
 597
 598 void
 599 lambda_matrix_project_to_null (lambda_matrix B, int rowsize,
 600                                int colsize, int k, lambda_vector x)
 601 {
 602   lambda_matrix M1, M2, M3, I;
 603   int determinant;
 604
 605   /* Compute c(I-B^T inv(B B^T) B) e sub k.  */
 606
 607   /* M1 is the transpose of B.  */
 608   M1 = lambda_matrix_new (colsize, colsize);
 609   lambda_matrix_transpose (B, M1, rowsize, colsize);
 610
 611   /* M2 = B * B^T */
 612   M2 = lambda_matrix_new (colsize, colsize);
 613   lambda_matrix_mult (B, M1, M2, rowsize, colsize, rowsize);
 614
 615   /* M3 = inv(M2) */
 616   M3 = lambda_matrix_new (colsize, colsize);
 617   determinant = lambda_matrix_inverse (M2, M3, rowsize);
 618
 619   /* M2 = B^T (inv(B B^T)) */
 620   lambda_matrix_mult (M1, M3, M2, colsize, rowsize, rowsize);
 621
 622   /* M1 = B^T (inv(B B^T)) B */
 623   lambda_matrix_mult (M2, B, M1, colsize, rowsize, colsize);
 624   lambda_matrix_negate (M1, M1, colsize, colsize);
 625
 626   I = lambda_matrix_new (colsize, colsize);
 627   lambda_matrix_id (I, colsize);
 628
 629   lambda_matrix_add_mc (I, determinant, M1, 1, M2, colsize, colsize);
 630
 631   lambda_matrix_get_column (M2, colsize, k - 1, x);
 632
 633 }
 634
 635 /* Multiply a vector VEC by a matrix MAT.
 636    MAT is an M*N matrix, and VEC is a vector with length N.  The result
 637    is stored in DEST which must be a vector of length M.  */
 638
 639 void
 640 lambda_matrix_vector_mult (lambda_matrix matrix, int m, int n,
 641                            lambda_vector vec, lambda_vector dest)
 642 {
 643   int i, j;
 644
 645   lambda_vector_clear (dest, m);
 646   for (i = 0; i < m; i++)
 647     for (j = 0; j < n; j++)
 648       dest[i] += matrix[i][j] * vec[j];
 649 }
 650
 651 /* Print out an M x N matrix MAT to OUTFILE.  */
 652
 653 void
 654 print_lambda_matrix (FILE * outfile, lambda_matrix matrix, int m, int n)
 655 {
 656   int i;
 657
 658   for (i = 0; i < m; i++)
 659     print_lambda_vector (outfile, matrix[i], n);
 660   fprintf (outfile, "\n");
 661 }
 662