1 /* Integer matrix math routines
2 Copyright (C) 2003, 2004, 2005 Free Software Foundation, Inc.
3 Contributed by Daniel Berlin <dberlin@dberlin.org>.
5 This file is part of GCC.
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
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
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
23 #include "coretypes.h"
29 static void lambda_matrix_get_column (lambda_matrix
, int, int,
32 /* Allocate a matrix of M rows x N cols. */
35 lambda_matrix_new (int m
, int n
)
40 mat
= GGC_NEWVEC (lambda_vector
, m
);
42 for (i
= 0; i
< m
; i
++)
43 mat
[i
] = lambda_vector_new (n
);
48 /* Copy the elements of M x N matrix MAT1 to MAT2. */
51 lambda_matrix_copy (lambda_matrix mat1
, lambda_matrix mat2
,
56 for (i
= 0; i
< m
; i
++)
57 lambda_vector_copy (mat1
[i
], mat2
[i
], n
);
60 /* Store the N x N identity matrix in MAT. */
63 lambda_matrix_id (lambda_matrix mat
, int size
)
67 for (i
= 0; i
< size
; i
++)
68 for (j
= 0; j
< size
; j
++)
69 mat
[i
][j
] = (i
== j
) ? 1 : 0;
72 /* Return true if MAT is the identity matrix of SIZE */
75 lambda_matrix_id_p (lambda_matrix mat
, int size
)
78 for (i
= 0; i
< size
; i
++)
79 for (j
= 0; j
< size
; j
++)
95 /* Negate the elements of the M x N matrix MAT1 and store it in MAT2. */
98 lambda_matrix_negate (lambda_matrix mat1
, lambda_matrix mat2
, int m
, int n
)
102 for (i
= 0; i
< m
; i
++)
103 lambda_vector_negate (mat1
[i
], mat2
[i
], n
);
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. */
110 lambda_matrix_transpose (lambda_matrix mat1
, lambda_matrix mat2
, int m
, int n
)
114 for (i
= 0; i
< n
; i
++)
115 for (j
= 0; j
< m
; j
++)
116 mat2
[i
][j
] = mat1
[j
][i
];
120 /* Add two M x N matrices together: MAT3 = MAT1+MAT2. */
123 lambda_matrix_add (lambda_matrix mat1
, lambda_matrix mat2
,
124 lambda_matrix mat3
, int m
, int n
)
128 for (i
= 0; i
< m
; i
++)
129 lambda_vector_add (mat1
[i
], mat2
[i
], mat3
[i
], n
);
132 /* MAT3 = CONST1 * MAT1 + CONST2 * MAT2. All matrices are M x N. */
135 lambda_matrix_add_mc (lambda_matrix mat1
, int const1
,
136 lambda_matrix mat2
, int const2
,
137 lambda_matrix mat3
, int m
, int n
)
141 for (i
= 0; i
< m
; i
++)
142 lambda_vector_add_mc (mat1
[i
], const1
, mat2
[i
], const2
, mat3
[i
], n
);
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. */
150 lambda_matrix_mult (lambda_matrix mat1
, lambda_matrix mat2
,
151 lambda_matrix mat3
, int m
, int r
, int n
)
156 for (i
= 0; i
< m
; i
++)
158 for (j
= 0; j
< n
; j
++)
161 for (k
= 0; k
< r
; k
++)
162 mat3
[i
][j
] += mat1
[i
][k
] * mat2
[k
][j
];
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. */
171 lambda_matrix_get_column (lambda_matrix mat
, int n
, int col
,
176 for (i
= 0; i
< n
; i
++)
177 vec
[i
] = mat
[i
][col
];
180 /* Delete rows r1 to r2 (not including r2). */
183 lambda_matrix_delete_rows (lambda_matrix mat
, int rows
, int from
, int to
)
189 for (i
= to
; i
< rows
; i
++)
190 mat
[i
- dist
] = mat
[i
];
192 for (i
= rows
- dist
; i
< rows
; i
++)
196 /* Swap rows R1 and R2 in matrix MAT. */
199 lambda_matrix_row_exchange (lambda_matrix mat
, int r1
, int r2
)
208 /* Add a multiple of row R1 of matrix MAT with N columns to row R2:
209 R2 = R2 + CONST1 * R1. */
212 lambda_matrix_row_add (lambda_matrix mat
, int n
, int r1
, int r2
, int const1
)
219 for (i
= 0; i
< n
; i
++)
220 mat
[r2
][i
] += const1
* mat
[r1
][i
];
223 /* Negate row R1 of matrix MAT which has N columns. */
226 lambda_matrix_row_negate (lambda_matrix mat
, int n
, int r1
)
228 lambda_vector_negate (mat
[r1
], mat
[r1
], n
);
231 /* Multiply row R1 of matrix MAT with N columns by CONST1. */
234 lambda_matrix_row_mc (lambda_matrix mat
, int n
, int r1
, int const1
)
238 for (i
= 0; i
< n
; i
++)
239 mat
[r1
][i
] *= const1
;
242 /* Exchange COL1 and COL2 in matrix MAT. M is the number of rows. */
245 lambda_matrix_col_exchange (lambda_matrix mat
, int m
, int col1
, int col2
)
249 for (i
= 0; i
< m
; i
++)
252 mat
[i
][col1
] = mat
[i
][col2
];
257 /* Add a multiple of column C1 of matrix MAT with M rows to column C2:
258 C2 = C2 + CONST1 * C1. */
261 lambda_matrix_col_add (lambda_matrix mat
, int m
, int c1
, int c2
, int const1
)
268 for (i
= 0; i
< m
; i
++)
269 mat
[i
][c2
] += const1
* mat
[i
][c1
];
272 /* Negate column C1 of matrix MAT which has M rows. */
275 lambda_matrix_col_negate (lambda_matrix mat
, int m
, int c1
)
279 for (i
= 0; i
< m
; i
++)
283 /* Multiply column C1 of matrix MAT with M rows by CONST1. */
286 lambda_matrix_col_mc (lambda_matrix mat
, int m
, int c1
, int const1
)
290 for (i
= 0; i
< m
; i
++)
291 mat
[i
][c1
] *= const1
;
294 /* Compute the inverse of the N x N matrix MAT and store it in INV.
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.
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
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. */
309 static int lambda_matrix_inverse_hard (lambda_matrix
, lambda_matrix
, int);
312 lambda_matrix_inverse (lambda_matrix mat
, lambda_matrix inv
, int n
)
325 det
= (a
* d
- b
* c
);
337 return lambda_matrix_inverse_hard (mat
, inv
, n
);
340 /* If MAT is not a special case, invert it the hard way. */
343 lambda_matrix_inverse_hard (lambda_matrix mat
, lambda_matrix inv
, int n
)
350 temp
= lambda_matrix_new (n
, n
);
351 lambda_matrix_copy (mat
, temp
, n
, n
);
352 lambda_matrix_id (inv
, n
);
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
357 for (j
= 0; j
< n
; j
++)
361 /* Make every element in the current row positive. */
362 for (i
= j
; i
< n
; i
++)
365 lambda_matrix_col_negate (temp
, n
, i
);
366 lambda_matrix_col_negate (inv
, n
, i
);
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
)
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
);
377 for (i
= j
+ 1; i
< n
; i
++)
381 factor
= -1 * row
[i
];
385 lambda_matrix_col_add (temp
, n
, j
, i
, factor
);
386 lambda_matrix_col_add (inv
, n
, j
, i
, factor
);
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. */
396 for (j
= n
- 1; j
>= 0; j
--)
403 /* The matrix must not be singular. */
404 gcc_assert (diagonal
);
406 determinant
= determinant
* diagonal
;
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
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
);
418 row
[j
] = diagonal
= 1;
421 /* Sweep the lower triangle column wise. */
422 for (i
= j
- 1; i
>= 0; i
--)
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
);
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
442 lambda_matrix_hermite (lambda_matrix mat
, int n
,
443 lambda_matrix H
, lambda_matrix U
)
446 int i
, j
, factor
, minimum_col
;
448 lambda_matrix_copy (mat
, H
, n
, n
);
449 lambda_matrix_id (U
, n
);
451 for (j
= 0; j
< n
; j
++)
455 /* Make every element of H[j][j..n] positive. */
456 for (i
= j
; i
< n
; i
++)
460 lambda_matrix_col_negate (H
, n
, i
);
461 lambda_vector_negate (U
[i
], U
[i
], n
);
465 /* Stop when only the diagonal element is nonzero. */
466 while (lambda_vector_first_nz (row
, n
, j
+ 1) < n
)
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
);
472 for (i
= j
+ 1; i
< n
; i
++)
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
);
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".
486 Ref: Algorithm 2.1 page 33 in "Loop Transformations for
487 Restructuring Compilers" Utpal Banerjee. */
490 lambda_matrix_right_hermite (lambda_matrix A
, int m
, int n
,
491 lambda_matrix S
, lambda_matrix U
)
495 lambda_matrix_copy (A
, S
, m
, n
);
496 lambda_matrix_id (U
, m
);
498 for (j
= 0; j
< n
; j
++)
500 if (lambda_vector_first_nz (S
[j
], m
, i0
) < m
)
503 for (i
= m
- 1; i
>= i0
; i
--)
507 int sigma
, factor
, a
, b
;
511 sigma
= (a
* b
< 0) ? -1: 1;
514 factor
= sigma
* (a
/ b
);
516 lambda_matrix_row_add (S
, n
, i
, i
-1, -factor
);
517 lambda_matrix_row_exchange (S
, i
, i
-1);
519 lambda_matrix_row_add (U
, m
, i
, i
-1, -factor
);
520 lambda_matrix_row_exchange (U
, i
, i
-1);
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".
531 Ref: Algorithm 2.2 page 36 in "Loop Transformations for
532 Restructuring Compilers" Utpal Banerjee. */
535 lambda_matrix_left_hermite (lambda_matrix A
, int m
, int n
,
536 lambda_matrix S
, lambda_matrix V
)
540 lambda_matrix_copy (A
, S
, m
, n
);
541 lambda_matrix_id (V
, m
);
543 for (j
= 0; j
< n
; j
++)
545 if (lambda_vector_first_nz (S
[j
], m
, i0
) < m
)
548 for (i
= m
- 1; i
>= i0
; i
--)
552 int sigma
, factor
, a
, b
;
556 sigma
= (a
* b
< 0) ? -1: 1;
559 factor
= sigma
* (a
/ b
);
561 lambda_matrix_row_add (S
, n
, i
, i
-1, -factor
);
562 lambda_matrix_row_exchange (S
, i
, i
-1);
564 lambda_matrix_col_add (V
, m
, i
-1, i
, factor
);
565 lambda_matrix_col_exchange (V
, m
, i
, i
-1);
572 /* When it exists, return the first nonzero row in MAT after row
573 STARTROW. Otherwise return rowsize. */
576 lambda_matrix_first_nz_vec (lambda_matrix mat
, int rowsize
, int colsize
,
582 for (j
= startrow
; (j
< rowsize
) && !found
; j
++)
585 && (lambda_vector_first_nz (mat
[j
], colsize
, startrow
) < colsize
))
594 /* Calculate the projection of E sub k to the null space of B. */
597 lambda_matrix_project_to_null (lambda_matrix B
, int rowsize
,
598 int colsize
, int k
, lambda_vector x
)
600 lambda_matrix M1
, M2
, M3
, I
;
603 /* Compute c(I-B^T inv(B B^T) B) e sub k. */
605 /* M1 is the transpose of B. */
606 M1
= lambda_matrix_new (colsize
, colsize
);
607 lambda_matrix_transpose (B
, M1
, rowsize
, colsize
);
610 M2
= lambda_matrix_new (colsize
, colsize
);
611 lambda_matrix_mult (B
, M1
, M2
, rowsize
, colsize
, rowsize
);
614 M3
= lambda_matrix_new (colsize
, colsize
);
615 determinant
= lambda_matrix_inverse (M2
, M3
, rowsize
);
617 /* M2 = B^T (inv(B B^T)) */
618 lambda_matrix_mult (M1
, M3
, M2
, colsize
, rowsize
, rowsize
);
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
);
624 I
= lambda_matrix_new (colsize
, colsize
);
625 lambda_matrix_id (I
, colsize
);
627 lambda_matrix_add_mc (I
, determinant
, M1
, 1, M2
, colsize
, colsize
);
629 lambda_matrix_get_column (M2
, colsize
, k
- 1, x
);
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. */
638 lambda_matrix_vector_mult (lambda_matrix matrix
, int m
, int n
,
639 lambda_vector vec
, lambda_vector dest
)
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
];
649 /* Print out an M x N matrix MAT to OUTFILE. */
652 print_lambda_matrix (FILE * outfile
, lambda_matrix matrix
, int m
, int n
)
656 for (i
= 0; i
< m
; i
++)
657 print_lambda_vector (outfile
, matrix
[i
], n
);
658 fprintf (outfile
, "\n");